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should  be  such  that  the  quantities  involved  are  all  smoothly  varying  functions. 

Numerical  simulation  can  contribute  to  this  goal  in  two  ways; 

(1)  It  can  provide  "benchmark"  problems  against  which  approximations  can  be 
tested.  (The  attack  on  these  problems  can  be  regarded  as  a  precision 
experiment . )  This  report  discusses  particular  problems  of  this  kind. 

(2)  It  can  help  to  provide  theoretical  insight  into  the  nature  of  turbulence. 
(For  example,  numerical  experiments  are  possible  that  are  not  feasible  in 
the  laboratory.)  Hopefully,  universal  features  will  be  found  and  then 
approximate  methods  can  be  developed. 

For  three-dimensional  problems  «t  is  shown  that  the  Reynolds  number  of  prob¬ 
lems  that  can  be  treated,  even  by  the  new  generation  of  computers,  is  rather 
limited  (R  <  10  to  4  X  10  ).  ’’he  Illiac  IV  is  especially  restricted  in  treating 
turbulence  problems.  It  has  a  relatively  small  fast  memory.  For  the  problems  of 
interest,  there  are  many  variables  and  the  interactions  are  nonlocal— i . e. , 
quantities  to  be  computed  at  a  point  depend  on  ones  at  all  other  points.  This  does 
not  mean  that  important  problems  cannot  be  simulated.  Transition  from  laminar  to 
turbulent  flow  and  the  behavior  at  low  Reynolds  number  above  transition  are  of 
interest .  Some  of  these  can  be  treated  with  existing  computers  (the  IBM  360/195 
and  the  CDC  7600).  However,  it  appears  unlikely  that  a  major  breakthrough  will 
necessarily  result  from  the  calculations.  But  even  these  calculations  may  serve 
to  disprove  some  approximation  methods. 

In  the  somewhat  less  realisfc  case  of  two  dimensions,  the  situation  is 
slightly  different.  Calculations  already  done  show  a  tendency  toward  Reynolds- 
number  independence  of  macroscopic  features.  If  confirmed  by  the  results  of  the 
higher  Reynolds-number  calculations  (made  possible  with  the  new  computers)  this 
could  be  very  significant.  (The  caveat  about  the  Illiac  IV  still  applies.) 

Finally,  we  consider  what  computer  capability  is  needed  to  mateiially  change 
the  above  picture,  and  sketch  how  this  could  possibly  be  achieved  and  at  what  cost. 
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ABSTRACT 


The  question  considered  is  whether  now  is  an  appropriate  time  for 
a  direc^  numerical  simulation  of  turbulence,  with  particular  emphasis 
on  the  use  of  the  Illiac  IV. 

First,  it  should  be’ emphasized  that  the  use  of  numerical  simulation 

as  a  general  tool  for  turbulence  calculations  is  out  of  the  question. 

To  resolve  all  scales  of  motion  when  the  Reynolds  number  is  R  requires 

9/4 

following  on  the  order  of  R  independent  variables  in  time,  and  on 
3 

the  order  of  R  JfcnR  arithmetic  operations  are  needed.  Since  many  problems 

6 

of  interest  involve  R  >  10  ,  it  is  clear  that  there  are  situations  well 
beyond  the  capabilities  of  exis-ing  or  projected  computers. 

The  realistic  goal  of  practical  turbulence  theory  is  to  achieve 
approximate  equations  for  mean  quantities.  These  equations,  even  though 
complicated  looking,  should  be  such  that  the  quantities  involved  are  all 
smoothly  varying  functions. 

Numerical  simulation  can  contribute  to  this  goal  in  two  ways: 

(1)  It  can  provide  "benchjnark"  problems  against  which  approxi¬ 
mations  can  be  tested.  (The  attack  on  these  problems  can 

be  regarded  as  a  precision  experiment.)  This  report  discusses 
particular  problems  of  this  kind. 

(2)  It  can  help  to  provide  theoretical  insight  into  the  nature 
of  turbulence.  (For  example,  numerical  experiments  are 
possible  that  are  not  feasible  in  the  laboratory.)  Hope¬ 
fully,  universal  features  will  be  found  and  then  approximate 
methods  can  be  developed. 

For  three-dimensional  problems  it  is  shown  that  the  Reynolds  number 
of  problems  that  can  be  treated,  even  by  the  new  generation  of  computers, 


ii 


is  rather  limited  (R  <  104  to  4  X  104).  The  Illiac  IV  is  especially 
restricted  in  treating  turbulence  problems.  It  has  a  relatively  small 
fast  memory.  For  the  problems  of  interest,  there  are  many  variables 
and  the  interactions  are  nonlocal— i .e. ,  quantities  to  be  computed  at 
a  point  depend  on  ones  at  all  other  points.  This  does  not  mean  that 
important  problems  cannot  be  simulated.  Transition  from  laminar  to 
turbulent  flow  and  the  behavior  at  low  Reynolds  number  above  transition 
are  of  interest.  Some  of  these  can  be  treated  with  existing  computers 
(the  IBM  360/195  and  the  CDC  7600).  However,  it  appears  unlikely  that 
a  major  breakthrough  will  necessarily  result  from  the  calculations. 

But  even  these  calculations  may  serve  to  disprove  some  approximation 
methods. 

In  the  somewhat  less  realistic  case  of  two  dimensions,  the  situa¬ 
tion  is  slightly  different.  Calculations  already  done  show  a  tendency 
toward  Reynolds -number  independence  of  macroscopic  features.  If 
confirmed  by  the  results  of  the  higher  Reynolds-number  calculations 
(made  possible  with  the  new  computers)  this  could  be  very  significant. 

(The  caveat  about  the  Illiac  IV  still  applies.) 

Finally,  we  consider  what  computer  capability  is  needed  to  materially 
change  the  above  picture,  and  sketch  how  this  could  possibly  be  achieved 
and  at  what  cost. 
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I  INTRODUCTION 


Turbulent  flows  play  a  major  role  in  many  problems  of  great 
practical  interest.  Among  these  are  flow  in  pipes  (industrial  appli¬ 
cations),  wakes  (reentry  vehicles,  submarines),  boundary  layers 
(separation,  transition),  meteorology,  airplane  designs,  and  nuclear 
explosions.  Unfortunately,  even  at  this  date  our  understanding  of 
turbulence  is  severely  limited.  Qualitative  features  are  known  as 
are  some  general  predictions  (based  usually  on  some  reasonable 
assumptions  and  some  dimensional  analysis) .  Practical  problems  are 
usually  treated  by  phenomenological  equations.  Tne  range  of  validity 
of  these  is  either  unknown  or  known  to  be  restricted  to  situations 
close  to  those  in  which  the  phenomenological  constants  are  determined. 

Here  the  question  is  considered  as  together  the  advent  of  a  new 
generation  of  computers  (Illiac  IV,  Star,  ASC)  makes  feasible  a  numerical 
attack  on  turbulence  by  direct  integration  of  the  Navier-Stokes  equations 

It  is  concluded  that,  while  some  useful  problems  can  (and  should) 
be  done,  it  does  not  seem  that  a  fundamental  breakthrough  is  perceivable. 
This  statement  is  particularly  true  for  the  Illiac  IV  because  of  memory 
limitations . 

This  should  be  understood  in  the  following  context: 

(1)  Presumably  most  practical  problems  are  so  complex,  both 
geometrically  and  physically,  that  we  do  not  expect  a 
direct  solution  to  le  possible  from  first  principles. 

v2)  The  main  obstacle  t >  numerical  solution  of  the  Navier-Stokes 
equation  lies  in  the  fact  that  when  there  is  turb»lence  we 
have  interesting  and  important  motions  on  many  different 
scale-  (both  space  and  time).  Tn  particular,  the  larger 
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the  Reynolds  number  the  smaller  are  the  scales  that  need  be 
considered.  Clearly,  for  any  computer  there  will  be  some 
Reynolds  number  beyond  which  the  computation  of  solutions 
is  not  reasonable.  We  do  not  expect  to  find  asymptotic 
behaviors  numerically,  but  rather  hope  that  the  numerics 
will  le  a  guide  to  an  asymptotic  theory. 


Accordingly,  we  phrase  the  question  as  follows:  Are  there  problems 
of  interest  that  can  now  be  treated  that  were  previously  inaccessible  to 
direct  numerical  solution?  By  "interest"  we  mean  either: 


They  are  benchmark"  problems:  The  results  should  be  veri¬ 
fiable  by  ’cperiments .  Assuming  agreement  is  achieved,  one 
then  wouiu  have  a  much  more  detailed  description  than  can 
ever  be  found  from  experiment.  (For  example,  one  could  cal¬ 
culate  as  many  correlation  functions  and  their  transforms 
as  desired.)  These  would  then  be  tests  to  which  approximation 
methods  would  be  subject. 

The  results  give  theoretical  insight.  For  example,  is  the 
transition  from  laminar  to  turbulent  flow  as  Reynolds  number 
increases  due  to  successively  more  instabilities  arising? 

Is  our  picture  of  energy  and  vorticity  cascade  correct? 
Alternatively,  a  better  understanding  of  these  effects  can 
be  obtained  by  varying  boundary  conditions  in  either  wave 
number  or  physical  sptce.  Here  the  calculations  are 

which  are  either  impossible  or  very  expensive 
to  do  in  a  laboratory. 


Essentially  what  we  would  like  to  do  is 


Determine  whether  there  are  indeed  universal  features  of 
turbulence . 

If  there  are  such,  can  we  develop  equations  which  describe 
the  average  flows  without  requiring  the  resolution  of  the 
fine  scales  of  the  turbulent  motion? 


To  get  some  id  x  of  the  barriers  encountered,  we  note  that  for  a 
problem  in  d  dimensions  with  N  degrees  of  freedom  in  each  linear 
dimension,  the  number  of  arithmetic  operations  to  be  done  is  on  the 
order  of 


where  a  and  b  ~  103.  As  a  rule  of  thumb,  we  need  N  ~  R^^/10  where  R 
is  the  Reynolds  number .  Two  things  stand  out :  For  a  given  computer, 
much  more  can  be  done  in  two  dimensions  than  three.  To  significantly 
increase  the  Reynolds  numbers  for  which  one  can  calculate  (say  by  a 
factor  of  two  to  ten)  requires  a  huge  increase  in  computer  capability. 
(This  is  particularly  true  in  three  dimensions.) 

However,  it  does  seem  that  some  interesting  problems  might  now  be 
done.  The  estimates  made  below  are  based  primarily  on  the  characteristics 
of  the  I  Iliac  IV  and  are  admittedly  somewhat  optimistic.  While  each  com- 
P’.ter  has  its  own  limitations  and  advantages  it  is  not  thought  that  the 
situation  would  be  radically  different  for  other  computers  of  the  new 
generation.  (As  mentioned  above,  the  Illiac  IV  does  have  memory  problems 
that  may  prove  intractable.) 

Problems  that  might  be  approached  are: 

(1)  The  behavior  of  flows  as  the  Reynolds  number  is  gradually 
increased  above  the  critical  value  in  Poisseuille  flow  in 
channels  and  pipes,  and  Couette  flow  in  channels  and 
between  rotating  cylinders. 

(2)  Turbulent  flows  in  pipes  and  channels  for  R  ~  I04. 

(There  are  many  experimental  results  available  for  these  two  classes 
or  problems.) 

(3)  Problems  of  homogeneous  isotropic  turbulence.  These  would 
be  primarily  numerical  expc  riments — to  enlarge  theoretical 
understanding. 

(4)  Much  higher  Reynolds  numbers  in  two  dimensions  can  be 
investigated.  This  is  of  interest  theoretically  and  has 
practical  use  for  atmospheric  problems. 

(5)  Problems  involving  turbi'lent  flows  in  boundary  layers— -say 
over  a  flat  plate.  This  is  something  that  probably  should 
be  one  of  the  last  to  be  considered.  The  difficulty  is  that 
to  treat  these  numerically  one  needs  to  put  in  some  (perhaps 
very  reasonable)  assumptions  as  to  the  behavior  of  the  flows 
at  large  distance  from  the  flow.  The  x*esults  might  merely 
reflect  these  assumptions  and  not  the  correct  physics. 


Some  important  questions  will  not  be  resolved  by  such  calculations. 

-5/3  ,  . 

The  existence  of  an  inertial  range  and  Kolmogoroff  k  law  (or  pro¬ 

posed  modifications)  would  not  be  seen  in  three  dimensions.  (However, 
the  analogous  k  might  be  observed  for  two-dimensional  problems.)  The 
interesting  questions  of  intermittency  would  also  not  be  particularly 
clarified. 

What  can  one  hope  to  be  the  conclusions  from  such  calculations? 
First,  we  might  verify  the  expectation  that  on  small  scales  turbulence 
is  essentially  universal— i  .e  . ,  rather  independent  of  the  large  scale 
structure.  Second,  one  might  check  previous  calculations  which  suggest 
that  at  sufficiently  high  Reynolds  numbers  the  macroscopic  features  of 
a  flow  are  rather  Reynolds-number  independent.  From  numerical  experi¬ 
ments,  new  methods  for  treating  sub-grid  scale  features  of  flow  may  be 
suggested.  Thes  (and  existing)  methods  can  then  be  tested  by  the 
calculations. 

It  may  be  noted  that  compressible  problems  have  not  been  mentioned. 
This  is,  of  course,  not  because  they  are  uninteresting.  Indeed,  some 
methods  of  calculation  exploit  compressibility.  However,  there  is  even 
less  understanding  of  turbulence  in  the  compressible  case.  In  light 
of  our  expectation  that  a  numerical  program  is,  at  least  in  part,  to 
extend  theoretical  understanding,  it  seems  best  to  start  where  we  feel 
we  now  have  the  most  knowledge. 

If  a  numerical  program  is  undertaken  (say  on  I Iliac  IV)  it  should 
go  forward  to  begin  with  at  a  rather  moderate  pace.  First  an  attempt 
should  be  made  to  see  how  a  program  might  be  written  for  some  simple 
problem.  It  seems  clear  from  our  estimates  that  problems  of  interest 
will  strain  the  capability  of  the  Illiac  as  it  is  expected  to  be  in 
the  next  year  or  so.  While  we  have  tried  to  anticipate  problems  that 
might  arise  due  to  pecularities  of  the  Illiac,  there  may  be  additional 


ones  which  rear  their  heads  when  a  detailed  program  is  attempted.  In 
particular  it  seems  very  likely  that  the  computation  scheme  that  should 
be  chosen  may  well  be  determined  by  properties  of  the  machine.  For 
example,  input-output  times  and  memory  capacity  may  well  dictate  a 
scheme  that  requires  somewhat  more  arithmetic  operations  than  a  perhaps 
otherwise  optimum  method. 

After  a  program  has  been  written,  numerous  tests  would  need  to  be 
performed.  These  might  be  other  (smaller)  calculations  done  in  two 
dimensions  and  also  on  homogeneous  isotropic  turbulence  in  three  dimen¬ 
sions.  The  standard  color-cone  problem  and  the  Taylor-Green  problem 
may  also  be  used. 

It  is  also  thought  that  the  whole  program  should  not  be  allocated 
to  a  single  contractor,  and  there  should  be  an  outside  advisory  group. 

The  reasoning  here  is  that  since  the  hope  is  to  gain  insight  as  the 
program  progresses,  one  would  expect  the  selection  of  succeeding 
problems  to  be  determined  by  what  has  baen  already  fcund. 

In  the  following  sections  we  amplify  and  sketch  the  reasoning  for 
the  above  remarks.  Section  II  briefly  describes  the  present  status  of 
turbulence  theory  and  experiment.  Section  III  describes  calculations 
which  have  been  done  and  those  that  might  be  done.  Some  general 
features  of  computation  procedures  for  numerical  simulation  of  turbu¬ 
lence  are  described  in  Section  IV.  Specific  details  of  the  Illiac  IV, 
which  determine  what  calculations  might  be  done  and  how  the  choice  of 
procedures  may  be  restricted,  are  discussed  in  Section  V.  Some  estimates 
of  computing  time  are  given  there.  Problems  of  interest  that  could  be 
done,  and  how  they  might  form  a  program  of  work  on  the  Illiac,  are  given 
in  St  ction  VI . 

A  brief  consideration  on  the  design  of  a  computer  that  might  be 
specific  for  turbulence  problems  is  given  in  Appendix  A. 


II  TURBULENCE 


A.  Status  of  Theory 

Most  flows  in  nature  are  turbulent— i  .e. ,  even  if  they  start  with 
a  relatively  smooth  laminar  profile  they  rapidly  develop  a  random  and 
fluctuating  appearance.  Qualitative  understanding  is  readily  obtained 
from  the  Nevier-Stokes  equations  for  incompressible  flow: 


7^  v  =  -VP  + 


2 

vV  v 


(1) 


(2) 


where  v(r, t)  is  the  velocity  field,  P(r,t)  is  the  pressure  and  v  is 
the  kinematic  velocity.  Typical  values  at  STP  are  (in  square  centi¬ 
meters  per  second)  v  J  =  0.145,  v  =  0.011,  v  =  0.00116. 

air  H20  Hg 

Let  us  estimate  the  ratios  of  the  non-linear  terms  to  the  viscous 
terms--! *  e • , 

v  •  Vv 

R  ~  - 

vV  v 


If  V  is  a  typical  value  of  the  velocity  and  L  a  measure  of  the  scale 
over  which  the  variations  are  occurring,  we  see 


v 


In  most  cases  this  is  a  large  number.  For  example,  for  water  flowing 

4 

in  a  pipe  of  radius  10  cm  with  velocity  10  cm/s,  R  =-  10  .  For  air 


flowing  at  10  m/hr  over  a  structure  of  radius  100  meters,  R  ^  4  X  10  . 
Under  these  circumstances  one  might  assume  that  the  viscous  terms  are 
completely  negligible.  However,  let  us  note  that  tne  convective  term 
(v  •  Vv)  is  nonlirear.  Expressed  in  wave  number  space,  this  shows  that 
the  various  Fourier  components  are  coupled  together.  Even  if  we  start 
with  only  small  wave  numbers  the  coupling  will  bring  in  higher  wave 
numbers.  For  these  wave  numbers  the  viscous  term  is  larger.  Eventually 
it  will  even  dominate  the  convective  term,  and  convert  the  flow  energy 
to  heat.  Thus  we  have  the  following  picture  for  flows  at  large  Reynolds 
number:  Energy  inserted  at  small  wave  numbers  cascades  to  higher  wave 
numbers  due  to  the  nonlinear  terms.  Eventually  viscosity  dominates  and 
the  energy  is  dissipated.  The  random  fluctuating  behavior  seen  at  large 
Reynolds  numbers  is  then  due  to  the  fact  that  high  wave  numbers  (i.e., 
high  Fourier  components)  are  excited. 

An  interesting  description  of  the  onset  of  turbulence  has  been 
given  by  Landau:1*  Below  a  certain  critical  Reynolds  number  the  laminar 
flow  is  stable.  Slightly  above  the  critical  Reynolds,  any  small  per¬ 
turbation  will  grow  slowly  until  the  solution  is  the  original  laminar 
flow  plus  a  small,  sinusoidally  varying  contribution  with  phase  deter¬ 
mined  by  the  perturbation.  For  larger  R  the  separation  into  laminar 
and  sinusoidal  becomes  no  longer  meaningful.  We  then  have  a  periodic 
flow,  but  not  merely  simply  periodic.  As  R  increases  further  we  get 
a  new  instability  and  eventually  a  doubly  periodic  flow  with  two  arbi¬ 
trary  phases  (essentially  determined  by  the  usually  undescribable  per¬ 
turbation),  Proceeding  to  higher  and  higher  R  we  get  successively  more 
and  more  nultiple  periodic  motions  with  more  and  more  arbitrary  phases. 
Hence  the  randomness  of  the  motion. 

* 

References  are  listed  at  the  end  of  this  report. 


Now  fcr  practical  purposes  v  e  are  not  interested  in  all  the  details 


of  the  fluctuating  motion,  but  rather  in  certain  average  values.  (For 
the  experimentalist  these  tend  to  be  either  spatia^  or  temporal  averages. 
The  theorist,  whose  point  of  view  we  take,  thinks  in  terms  of  ensemble 
averages,  and  implicitly  an  ergodic  theorem.)  With  bars  denoting  the 
appropriate  averages  it  is  conventional  to  divide  the  flov'  into  its 
average  and  fluctuating  parts — i.e.,  v  =  v  +  v/.  Taking  averages  of 
Eqs.  (1)  and  (2)  gives 


6t 


v  +  v 

i  J 


v 

i 


(3) 


These  are  perfectly  nice  equations  for  v,  but  involve,  however,  the 

Thus,  we  have  more  unknowns  than  equations. 

Of  course,  we  can  obtain  equations  for  the  Reynolds  stress  by  multi¬ 
plying  the  original  Eqs.  (1)  and  (2)  by  v  and  then  averaging.  But  then 


Reynolds  stress,  v^v^ 


quantities  like  p^v*  and  v * v *v^  occur.  Again  more  unknowns  than  equa¬ 
tions.  Proceeding  to  higher  order  we  find  the  came  situation  at  each 
step.  Breaking  this  chain  is  the  closure  problem.  Almost  all  approxi¬ 
mation  methods  are  based  on  some  assumption  as  to  how  this  can  be  done. 
For  example,  the  way  the  Reynolds  stress  occurs  in  Eq .  (3)  suggests 
replacing  that  term  by  one  like  that  duo  to  the  kinematic  viscosity 
but  with  a  fictitious  effective  viscosity.  This  eddy  viscosity  could, 
of  course,  depend  on  r  and  v,  and  even  more  complicated  quantities  to 
be  discussed  later.  Alternatively,  assumptions  have  been  made  about 
higher  order  correlation  coefficients.  It  may  be  noted  that  the 
resulting  equations  are  then  frequently  more  complicated  looking  than 
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Eqs.  (1)  and  (2).  However,  they  are  hopefully  for  smooth,  y  varying 
functions  and  therefore  easier  to  use  in  calculatioi  s. 

What  can  we  say  about  the  Reynolds  stress  tensor  using  our 
qualitative  ideas  and  some  dimensxonal  analysis?  If  the  ideas  about 
cascading  in  wave  number  space  are  right  we  expect  that 

v'(r, t)v'(r+s, t) 

is  rather  independent  of  r  for  small  s.  Let  cp^(k)  bo  the  Fourier 
transform  with  respect  to  s-~i.e,, 

00 

-3  rt  rt  -lk*8  '  f  ^ 3 

(k)  =  ( 2 tt )  -  ~  v^(r,t)  V^(r+a,t)  d  s  . 

— <30 

The  trace  of  this  quantity  can  be  interpreted  as  the  kinetic  energy 
at  a  given  wave  number «  Indeed, 

v'(r,t)  v  |(r,  t)  =  3u2  *  JJJ  cp^(k)d3k 

Def ining 

E(k)  «  {  JT  qp  (k)  k2  dQ 

J J  ii  -k 

(the  ir  tegral  Is  over  a  spherical  shell  of  radius  k),  wo  havo 

CD 

2 

E(k)  dk  ■  3/2  u  ■  kinotic  energy  per  unit  mass. 

0 

Consider  the  following  simple  model:  We  have  s  stationary  state 

2 

with  turbulent  energy,  u  ,  being  introduced  over  a  macroscopic  length. 
L.  Dy  the  cascade  process,  this  is  then  transferred  to  higher  wave 
numbers.  A  measure  of  the  time  for  this  to  take  place  is  i/u.  Hence 


0 


the  rate  at  which  energy  is  fed  into  the  cascade  is  u  /(jfc/u)  =  u  /l. 

At  sufficiently  high  k  this  energy  is  dissipated  by  viscosity  at  a 

3 

rate  e  which  (if  we  are  to  have  stationarity )  must  be  e  =  u  /l.  In  the 
dissipation  region  the  only  relevant  parameters  are  presumably  the 
kinematic  viscosity  v,  and  the  dissipation  rate  e.  The  appropriate 
time  and  space  scales  relevant  are  then  those  which  can  be  formed 
from  these  quantities.  They  ure  the  Kolmogorof f 2  microscales  of  time 

1.  3  X 

t  =  (v/e)  ,  distance  k  =  (v  /eT  .  From  these  one  can  form  the  velocity 
X 

scale  v  a  (ve)4.  Now  let  us  look  at  ranges  of  k  which  are  considerably 
larger  than  1/t.  One  expects  E(k)  to  be  a  function  only  of  e,  l,  and  k. 
Then  on  dimensional  grounds 


2 

E(k)  a  u  t  F(kZ,kn)  (4a) 

or  alternatively 

2 

E(k)  a  v  x  C(kZ,kK)  .  (4b) 

Note:  The  various  time,  distance  and  velocity  scales  are  all  related 

to  the  Reynolds  number.  Thus 


u 


(5) 


Now  let  us  suppose  that  we  have  a  region  of  k  such  that: 

kt  »  1,  and  simultaneously  kx  «  1  .  (0) 

From  Kq ,  (4a)  wo  then  expect 

Kfc)  ~  u 2  £  HkZ.O) 


.0 


and  i*om  Eq  .  (4b) 


2 

E(k)  **  v  k  G(oo,kH) 

2 

Equating  this  then  shows  v  k  G(oe,kK)  is  independent  of  v.  This  can 
onlv  to  1 1  ue  if  G(co,kn)  =  a  (kn)  "  with  a  some  universal  constant. 

Therefore  we  expect  that  there  may  be  a  region  of  k  such  that8 

E(k)  =  o  k'5  3  .  f 


Remarkably  enough  there  are  geophysical  fl„ws  in  which  a  region  of 
E(k)  ol  this  form  is  not  inconsistent  with  the  data  (with  a  ~  ],*>), 

\'e  note  though,  that  this  requires  there  be  a  region  of  k  far  from 
l>oth  those  k  directly  affected  by  macroscopic  variations  and  much 
less  than  the  k  region  subject  to  significant  dissipation.  Further, 
the  inequalities  of  Eq .  (6)  must  hold.  Combining  these  we  note  a 
necessary  condition  for  '  ich  a  region  to  exist  is 

L  2 

k*  «  1  <N  k Jt  or  -  ~  R  ‘  »»  l 

K 

That  is,  the  Reynolds  number  must  be  very  large. 

Sinlar  dimensional  arguments3  urc  extremely  useful  in  obtaining 
estimates  of  the  velocity  profiles  in  the  vicinity  of  walls  (wc  will 
need  a  knowledge  of  this  behavior  in  the  following.  Consider  the 
Idealised  case  of  a  steady  plane-parallel  flow  of  fluid  moving  in  «ie 
x  direction  In  the  spaco  between  moving  rigid  walls  at  r.  ■  0  and  t  ■  h 
in  the  absence  of  a  mean  pressure  gradient.  The  x -component  of  Eq .  (21 
gives 
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T(Z)  = 


2» 

d_v^ 

2 

dz 


—  v  v  =  0 
dz  x  z 


V 


dz 


v  v  =  constant  =  t 
x  z  o 


The  average  characteristics  of  the  flow  can  depend  only  on  the  shear 

stress  t  ,  kinematic  viscosity,  and  the  distances  z,  h-z.  It  is 
o 

conventional  to  introduce 

* 

u  =  /r 

o 


Then  dimensionally  we  expect 


V  (z)  =s  U  f 
X 


(h-z )u 


;)  ■ 


Now,  very  much  as  In  the  derivation  of  the  Kolmogoroff  spectrum,  we 
expect  that  for  a  sufficiently  la**ge  R  there  will  xist  a  region  of 
such  that 


zu 

—  finite, 


(h-z)u 

111  very  large 

V 


Then  approximately 


v  (z) 
x 


At  the  wall  v / 


0,  Then  Eq.  (8)  tolls  us 


dv 


dr. 


r«0 
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(8) 


Then,  as  the  first  term  in  a  Taylor  series  expansion,  we  have 


i  .  e  . , 


f(y)  =  y 


It  may  bo  noted  that  using  the  other  N-S  equations,  one  can  rhow  that 

4 

the  next  term  in  the  expansion  is  y  ,  Presumably  our  approximation  is 

♦  *  v 

good  out  to  some  z  of  the  order  of  z  =  —jr  I  i.e.,  the  linear  law  holds 
v  u 

out  to  a  z  =  a  -r  ,  From  experiment  one  finds  a  -  5  (this  region 

v 

0  £  z  £  a  -p-  is  called  the  viscous  sublayer!. 

When  we  are  far  out  of  the  viscous  sublayer  but  still  such  that 
u 

(h-z)  is  effectively  ,  we  can  also  get  an  expressior  for  f.  In 

such  regions  the  viscous  terms  are  negligible  compared  to  the  Reynolds 
dvx 

stresses.  Then  — —  can  only  depend  on  z  and  t  .  Therefore 


dz  <  z 

where  <  (the  von  Karman  constant)  is  presumably  universal.  Integration 
yields 


v  ( z )  *  — 

X  ' 


In 
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where  3  Is  a  new  universal  integration  constant.  Again  f  ror  experiment 
one  finds  *  O.i.  g  -  0.1.  What  are  the  limits  of  validity?  Clearly 
there  is  a  lower  limit  z  «  (const)  ./u*.  Kmpl rlcal ly ,  the  constant  is 
of  order  30.  We  then  have  the  following  picture.  For  z  ~  5v/u*  we  have 
a  linear  law.  lor  z  30  v/u  there  is  a  logarithmic  law.  In  between 
there  is  some  buffer  region  where  the  two  rur  match  up.  The  upper 


limit  of  the  logarithm  law  clearly  occurs  when  effects  of  the  other 
boundary  become  important . 


A  principal  object  of  theoretical  research  in  turbulence  for  the 
last  few  decades  has  been  isotropic  homogeneous  turbulence.  There  is 
a  basic  assumption  that  essential  features  of  turbulence  are  universal. 
This  is  in  agreement  with  our  above  qualitative  description  of  the 
cascade  process:  as  energy  goes  to  higher  wave  numbers,  memory  of  the 
initial  source  is  lost.  Further,  since  only  ensemble  averages  are 
eventually  desired,  a  strict  probabilistic  description  is  used.  For 
simplicity,  the  assumption  is  made  that  the  problem  is  homogeneous  and 
(usually)  isotropic.  Mathematically  the  formulation  is:4  Given  an 
infinite  fluid  described  by  Eqs.  (1)  and  (2),  and  given  the  initial 
value  of  the  velocity  as  a  random  function  of  position  (described  by 
certain  probability  laws  which  are  independent  of  position  and  direction), 
we  are  to  determine  the  probability  laws  that  describe  th<*  subsequent 
motion.  An  article  of  faith  is  included  here.  It  is  hoped  (conjectured) 
that  nonlinear  systems  with  a  large  number  of  degrees  of  freedom  will 
approach  a  statistical  state  which  is  essentially  independent  of  the 
initial  conditions. 

Clearly  there  is  no  exact  realization  of  isotropic  homogeneous 
turbulence  In  nature.  Indeed  some  features  realizable  In  shear  flows 
zre  certainly  absent  In  the  model.  However,  there  are  experiments  in 
wind  tunnels  (see  below)  In  which  the  turbulence  generated  doe®  appear 
to  be  approximately  homogeneous  and  Isotropic.  From  our  point  of  view, 
homogeneous  Isotropic  turbulence  should  probably  be  considered  as  an 
Interesting  model  from  which  one  hopes  to  get  theoretical  Insight. 

D  .  Status  of  Experiments 

Experimental  studies  of  turbulence  (Intensities,  spectra,  correla¬ 
tions,  etc.)  have  been  carried  out  for  about  the  last  50  year"  with 


electronic  instrumentation*  Although  many  of  tht  ie  experiments  were 
designed  to  study  flows  with  complicated  geometries  that  cannot  be 
readily  simulated  numerically,  some  of  them  provide  information  that 
can  be  used  to  check  the  validity  of  the  simulations.  Unfortunately 
most  experimentalists  have  attempted  to  obtain  as  large  Reynolds 
number  as  possible  in  order  to  study  the  asymptotic  state  of  the  flow, 
while  simulations  appear  to  be  possible,  at  present,  only  at  low  and 
moderate  Reynolds  numbers.  Thus  only  a  small  subset  of  the  existing 
turbulence  data  is  useful  for  comparison  with  simulation  results. 

The  measurements  referred  to  below  are  almost  exclusively  Eulerian. 
Due  to  the  obvious  experimental  difficulties,  almost  no  Lagrangian 
measurements  have  been  made  except  for  the  simplest  diffusion  experiments. 
One  of  the  advantages  of  simulations  is  the  ease  with  which  Lagrangian 
quantities — multiple  particle  correlation  coefficients,  for  example- 
can  be  calculated. 

1 .  Homogeneous — Isotropic  Turbulence 

The  simplest  flow  to  simulate— unbounded  homogeneous,  isotropic 
turbulence— is ,  unfortunately,  impossible  to  realize  exactly  in  an  ex¬ 
periment.  There  have  been  many  attempts  to  approximate  isotropic 
turbulence  experimentally.  Usually  these  involve  using  grids  in  wind 
tunnels ,  although  measurements  have  been  made  on  the  axis  of  a  circular 
free  Jet  and  measurements  of  the  small-scale  structure  have  been  made 
in  tidal  channels  and  in  the  atmospheric  boundary  layer.  The  geophy¬ 
sical  flows  do  not  produce  data  suitable  for  comparison  with  simulations 

because  R  is  large  (103  to  10^)  in  these  flc.s.  'Note:  R  ~  2R,2.) 

X 

The  wind  tunnel  and  Jet  flows  provide  more  useful  data,  for 

comparison  purposes,  with  R  in  the  range  of  10  to  200.  The  measure- 

K 

ments  generally  consist  of  turbulence  intensity,  spectra,  and  second 
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and  third  order  velocity  correlations.  Most  of  the  measurements  made, 
up  to  about  1961,  are  summarized  by  Batchelor4 and  Hinze.6 

In  recent  years,  furt  experiments  have  been  performed  in 

which  great  care  was  taken  to  make  the  flow  as  isotropic  as  possible. 

A  few  examples  of  these  are  described  below. 

(1)  Champagne,  Harris,  and  Corrsin’  carried  out  an  experi¬ 
ment  on  grid-generated  turbulence  in  a  wind  tunnel  at 
R\  =  130.  They  measured  turbulence  intensities, 

Reynolds  stress,  the  spectra  of  each  of  the  velocity 
components,  the  Taylor  microscale,  the  cross-spectrum 
of  u  and  v  and  various  second  order  space-time  correc¬ 
tions  of  u. 

(2)  Wyngard  and  Tenekes7  measured  the  probability  density, 
skewness  and  kurtosis  of  (bv'/bt)  and  (b2v  Vbt2)  as  well 
as  the  spectrum  of  (b2v'/bt2)  at  R^  =  200.  This  was 

in  a  mixing  layer. 

(3)  Comte-Bellot  and  Corrsinp  used  grid-generated  turbulence 
and  measured  turbulent  energy,  the  dissipation  rate  and 
second-order  correlation  functions  of  velocitv  at 

R^  =  35  to  70. 

(4)  Kuo  and  Corrsin,9  again  using  grid  turbulence,  measured 
the  probability  density  of  v',  (bv'/bt),  (b2v'/bt2),  u2, 
(bv'/bt)2  and  (b2vVbt2)  ai:l  the  kurtosis  of  (bv'/bt) 
and  (5  v  /frt  ) .  These  measurements  were  made  at  =  72 
and  830. 

These  measurements,  as  well  as  other  recent  work  and  the 
older  studies,  provide  a  wealth  of  experimental  information  on  "nearly 
isotropic  flow  which  can  be  compared  with  the  results  of  simulations. 

2 .  Pipe  and  Channel  Flow 

The  classic  experiments  are  those  of  Laufer10  at  R  =  50,000 
and  500,000  in  a  circular  pipe  and  in  a  wide  channel.  The  mean  profiles, 
turbulence  intensity,  the  Reynolds  stress,  turbulent  energy,  and  tur¬ 
bulent-energy-dissipation  rate  were  measured  at  various  distances  from 
the  wall  down  to  zu*/v  **  3.  The  energy  spectra  of  the  three  velocity 
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components  were  measurer!  at  several  distances  from  the  wall.  Seme 
third-order  correlations  at  a  poinc  and  some  second-order  spatial 
correlations  were  also  measured. 

Since  Laufer's  experiments,  other  workers  have  continued  to 
investigate  these  flows.  A  sample  of  some  of  these  follow: 

(1)  Bakewell  and  Lumley11  measured  the  mean  profile  for 
zu*/v  >  2;  the  energy  of  the  stream-wise  component  of 
the  turbulent  velocity,  the  spectra  and  the  second- 
order  correlation  of  u,  all  for  1  ^  zu*/v  ^  40.  These 
measurements  were  made  at  R  =  8700. 

(2)  Wallace,  Eckelmann,  and  Brodkey1"5  have  studied  channel 
flow  at  R  »  7150.  The  mean  profile  and  Reynolds  stress 
were  measured  in  the  range  1  ~  zu*/v  >  200,  as  well  as 
the  second-order  correlation  of  u  at  zu  /v  -  15. 

These  experiments  should  prove  the  most  valuable  for  comparison  with 
simulations  because  these  flows  can  be  simulated  ’exactly.” 

There  is  one  other  flow,  Couette  flow  between  rotating 
cylinders,  that  can  also  be  ’’exactly”  simulated.  However,  there  does 
not  seem  to  be  any  quantitative  study  of  turbulent  Couette  flow,  pos¬ 
sibly  because  of  the  difficulties  of  inserting  a  probe  that  does  not 
disturb  the  flow  too  much.  It  is  possible  that  laser-velocimeter 
techniques  could  be  applied  to  make  these  measurements. 

3 .  Boundary  Layers,  Wakes,  and  Jets 

There  is  an  enormous  literature  (see  Coles  and  Hirst,13  and 
Hinze6  and  Townsend  14  ior  summaries)  containing  measurements  of 
transition,  mean  profiles,  Reynolds  stresses,  correlation  functions, 
etc.,  for  boundary  layers,  wakes,  and  jets.  Because  of  the  difficult., 
in  simulating  these  flows,  these  experiments  are  probably  less  valuable 
for  comparison  purposes. 
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Ill  TURBTJD2NCE  CALCULATIONS 


A .  Present  Status  of  Turbulence  Simulations 

Turbulence  calculations,  using  closure  approximations  of  various 
types,  have  been  made  for  many  years.  By  suitable  adjustments  of  the 
t jrbulence  model  and  the  arbitrary  constants  contained  therein,  reason¬ 
able  agreement  with  experi'aent  has  been  obtained  for  particular  flows. 
None  of  these  models  appears  to  be  universal  and  it  is  not  even  known 
whether  or  not  a  universal  model  is  possible. 


A  few  direct  simulations  of  turbulence  have  been  carried  out  in 
recent  years.  These  were  carried  out  at  low  and  moderate  values  of  R 

A 

and  it  is  believed  that,  in  these  calculations,  the  resolution  was  fine 
enough  to  resolve  all  scales  down  to  the  Kolmogoroff  microscale.  The 
published  simulations  have  treated  homogeneous-isotropic  turbulence  in 
i'"o  or  three  dimensions. 


Lilly  15  and  Deem  and  Zabusky16  have  simulated  two-dimensional 
isotropic  turbulence  at  Reynolds  numbers  ranging  from  several  hundred 

to  a  few  thousand  (R  =  25  to  150).  Doubts  have  been  raised  about  the 

A 

adequacy  of  the  resolution  of  the  simulations  at  the  higher  values  of  R; 
however  there  is  no  doubt  that  the  resolution  was  adequate  for  R  ~  500 

(R.  £  60). 

A 

A  single  set  of  simulations  of  three-dimensional  isotropic  turbulence 
has  been  published  by  Orszag1^  and  Patterson, 18  using  both  finite  dif¬ 
ference  and  spectral  techniques.  In  these  calculations  R  =  20  to  40. 

A 

To  date  the  results  of  these  calculations  have  been  used  to  test 
various  theories  of  turbulence  in  both  two  and  three  dimensions.  These 


calculations  are  extremely  valuable  as  tests  of  theories  of  isotropic 
turbulence  because  homogeneous  isotropic  turbulence  cannot  be  exactly 
realized  in  the  laboratory  or  in  geophysical  flows.  In  particular,  the 
initial  conditions  of  the  simulations  can  be  varied  arbitrarily  and  any 
desired  information  about  the  flow  can  be  obtained  from  the  simulations 
with  relative  ease. 

Finally,  it  should  be  noted  that  a  number  of  workers  are  carrying 
out  simulations  that  are,  as  yet,  unpublished. 

B .  Calculations  of  Interest 

It  should  be  emphasized  that  a  direct  numerical  attack- -i  .  e . , 
starting  from  Eqs.  (1)  and  (2) — on  most  problems  of  practical  interest 
is  impossible  now  and  indeed  may  always  be  so.  The  reasons  are  the 
complicated  geometries  involved  and,  more  importantly,  the  range  of 
scales  that  must  be  described.  Below  we  indicate  that  in  three  uimen- 

sions  the  number  of  computations  to  be  done  for  one  time-dependent 

3  6  7 

problem  is  of  order  R  £nR.  Practically,  R  can  be  10  to  10  or  larger. 

The  goal  of  numerical  simulation  is  to  determine  universal  features  of 
turbulence.  If  such  are  found,  one  hopes  to  develop  approximate  equa¬ 
tions  for  the  average  flow  such  that  the  fine  sides  need  not  be  resolved. 
(Even  a  negative  result — that  there  are  some  features  of  different  flows 
that  are  not  universal — would  be  important.) 

Let  us  consider  then  what  problems  might  be  usefully  attacked 
numerically  with  the  new  generation  of  computers.  We  choose  two 
criteria . 

1 .  Benchmark  Problems 

These  p~e  problems  on  which  approximation  methods  could  be 
tested.  They  should  be  such  that  as  few  assumptions  as  possible  are 
put  into  the  calculation.  Existing  (or  possible)  experimental  results 


should  be  available.  Flows  in  pipes  and  channels  of  both  the  Poiseuille- 
and  Couette-type  are  of  this  nature.  So  is  Couette  flow  in  concentric 
cylinders.  Boundary-layer  flows  over  plates  are  less  useful  since  some 
extra  conjectures  as  to  the  behavior  at  infinity  are  necessary. 

2 .  Theoretical  Problems 

These  are  problems  whose  solution  may  lead  to  theoretical 
insight — and  then  to  practically  useful  approximations.  (Of  course, 
they  can  also  be  used  to  test  existing  approximations.)  Problems  of 
this  type  are: 

(1)  Studies  of  flows  just  above  the  critical  Reynolds  number: 
These  could  greatly  improve  our  understanding  of  the 
mechanics  of  transition.  An  example  is  the  transition 

in  Poiseuille  flow  in  pipes.  Linear-stability  theory19 
says  this  is  stable  for  all  R,  In  practice  the  transition 
occurs  at  R  ~  2200.  Clearly  the  Landau  picture1  of 
transition  discussed  above  cannot  be  correct  here.  The 
general  theory  of  nonlinear  instabilities  could  be  much 
clarified . 

(2)  Problems  of  isotropic  homogeneous  turbulence:  Presumably 
these  should  be  done  with  periodic  boundary  conditions. 
While  idealized  thesa  can  be  considered  as  numerical 
experiments.  For  e.  .ample,  if  one  has  results  a 

very  fine  grid,  one  could  see  how  well  this  could  be 
modeled  using  a  coarser  grid  with  approximate  boundary 
conditions  (in  wave  number  space). 

(3)  Two-dimensional  problems:  Though  also  somewhat  ideal, 
there  is  some  indication  that  the  atmospheric  problem 
is  of  this  type.  (There  is  some  disagreement  as  to 
whether  "turbulence"  can  be  two-dimensional  since  one 
method  of  vorticity  production--vortex  s tretching — is 
absent.  However,  adopting  the  view  that  turbulence 
refers  to  hydrodynamical  problems  that  require  a 
statistical  description,  one  concludes  that  there  is 

such  a  thing.)  Theoretically  the  two-dimensional  problems 
are  particularly  interesting  since,  as  we  shall  see, 
problems  with  much  more  interesting  Reynolds  numbers  can 
be  computed  than  are  possible  in  three  dimensions. 


Indeed,  there  are  indications  that  macroscopic  proper¬ 
ties  of  such  flows  are  insensitive  to  the  Reynolds 
number.  However,  these  simulations  may  not  have 
resolved  properly  all  necessary  scales  of  motion. 

For  accurate  simulation  the  enti  ^py-dissipat  ion 
spectrum  should  lie  well  within  the  resolvable  scales. 

(4)  Varying  boundary  conditions:  The  effects  of  boundaries 
can  be  studied  by  numerical  experiments.  Instead  of 
no-slip  conditions  in  pipes,  one  could  consider  other 
conditions  to  see  what  the  effects  are.  (These  are 
"experiments”  that  are  not  feasible  in  laboratories.) 

Of  course,  in  order  to  be  useful,  the  numerical  calculation 
must  be  done  for  problems  of  sufficiently  high  Rey.  'Ids  number.  For 
a  given  problem,  we  must  be  above  the  critical  value.  Also,  for  any 
material  increase  in  our  knowledge,  R  should  be  significantly  greater 
than  has  been  used  for  calculations  before. 

Some  very  interesting  questions  seem  beyond  the  range  of 
present  numerical  attack.  One  such  is  the  existence  of  an  equilibrium 
range  and  in  particular  an  inertial  subrange  where  the  Kolmogoroff 

spectrum  holds.  Theoretically  this  spectrum  is  a  little  peculiar. 

2/3 

We  see  that  Eq .  (7)  related  E(k),  (the  energy  density),  to  e  (the 
energy  dissipation  rate).  Since  e  is  itself  a  fluctuating  quantity, 
what  is  meant  here?  Presumably  an  average  is  meant  here.  But  E(k) 
is  also  an  average.  It  seems  strange  to  relate  one  average  to  a 
fractional  power  of  another.  Experimentally  there  is  also  some 
question,  Kolmogoroff ' s  assumption2  implies  not  only  the  5/3  spectrum 
but  also  that  other  statistical  quantities  of  the  small  scales  have 
universal  values  when  put  in  non-dimensional  form  with  e  and  v.  In 
particular,  in  the  inertial  range  there  should  be  universal  scaling 
with  e  alone.  Thus  the  skewness 


,  and  the  kurtosis 


Av  ;(r) |4 


where  Av;(r)  =  v;(x  +  r)  -  v'(x),  should  be  universal  constants.  How¬ 
ever,  measurements  suggest  that  skewness  and  kurtosis  rise  with  R. 

Attempts  to  accommodate  these  difficulties  give  rise  to  slight 
modifications  of  the  5/3  power  law.*'1'**  Also,  models**3  and  approximate 
calculations  tend  to  give  rise  to  other  laws  which  are  similar  to  this 
but  not  the  same. 


Can  we  resolve  the  question  by  numerical  simulation?  Probably 
today  we  cannot.  Thus,  for  the  existence  of  an  inertial  range,  we  saw 
that  one  must  have  a  range  of  k  such  that  simultaneously 


kn  «  1  and  kl  »  1  . 

Of  course  it  is  somewhat  arbitrary  as  to  what  these  inequalities  must 

be,  but  let  us  assume  100 :J  is  a  reasonable  value.  To  see  a  power  law 

one  would  need  several  decades  in  k.  Let  us  suppose  this  is  two.  At 

the  lower  end  (k  )  we  must  have 
o 


k  71  «=>  10 

o 


-4 


k  l 
o 


10 


and  at  the  upper  end  (k^)  »ve  have 


k  T|  10 


-2 


k  l  =  10 
1 


therefore 


1/ 7]  -  10 
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Indeed,  as  we  aoe  below,  Is  the  order  of  the  number  of  linear  grid 

G 

points  one  needs  for  a  calculation,  and  10  will  be  found  to  be  pro¬ 
hibitive. 

The  Reynolds  numbers  that  can  be  used  for  a  given  computer 

are  determined  by  the  resolution  r  luired.  In  the  case  of  Isotropic 

homogeneous  turbulence,  we  have  only  two  lengths— L,  a  periodicity 

length,  and  T),  the  Kolmogoroff  length  scale.  Somewhat  aoro  stringent 

conditions  are  act  In  contained  flows,  juch  as  pipes.  Here,  In  addition 

to  r  (the  pipe  radius)  and  T|,  we  have  the  thickness  of  the  viscous 
o 

sublayer  6  -  5  v/U  ,  and  the  thickness  of  the  region  over  which  there 

v 

it  a  change  from  the  linear  law  to  the  logarithmic  region,  4^  30  v/U  * 

• 

From  experiment  we  take  U  (the  friction  velocity)  and  U  the  turbulence 
U 

level  to  be  ~  with  U  the  moan  velocity  flow.  With  the  Reynolds 
20  a 

number  R  defined  by  R  ■  U  r  /v,  the  various  lengths  are  related  as: 

a  o 

10  100  500 

r  3/-1  *  v  R  '  b  R 

o  R 

Table  1  gives  these  rates  for  various  Interesting  Reynolds  numbers. 

Tablo  1 


RATES  TOR  SELECTED  REYNOLDS  NU101ERS 


R 

Tl/r 

o 

K/r 

b  o 

t  /r 
v  o 

2,000 

0.03 

0.25 

0.05 

5,000 

0.017 

0.10 

0.02 

10,000 

0.010 

0.05 

0.01 

20,000 

0.006 

0.025 

0.005 

40 , 000 

0.004 

0.125 

0.0025 

100,000 

0.002 

0.005 

0.001 

1,000,000 

0.0003 

0.0005 

0.0001 

From  this  table  sose  conclusion*  can  be  drawn.  For  ft  S  10*.  adequate 

resolution  la  Achieved  with  a  linear  grid  of  about  100  points; 

4  4  5 

ft  -  2  x  10  to  i  x  10  require*  about  200  points;  ft  -  10"  requires 

6 

about  500  to  1000  points,  *h lie  ft  -  10  requires  between  3000  and 
1 0 . OO0  |K>  Ints. 
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IV  COMPUTAT IOHAL  MLTHODS 


The  general  problem  la  lo  find  solutions  ol  l-qs.  (1)  and  ,i)  In 
«  region  V  with  appropriate  boundary  condition*  on  tho  bounding  surface 
S.  Many  scheme*  have  been  proposed  (and  used'  for  numerical  calcula¬ 
tion.  They  can  be  characterl red  in  many  different  ways. 

One  say  la  in  terms  of  Pierian  or  lagranglan  fora.  The  Lagranglan 
foraul  at  lona  are  characterised  by  a  coordinate  ayatca  that  aovea  with 
the  flo«.  Tor  f  lo»*  that  are  relatively  uncoapl  lea  it'd,  thl#  la  a  very 
attractive  method.  Tor  the  violent  f lows  of  a  turbulent  fluid  there 
are  aerloua  disadvantages  and  aoae  poaaible  advantage*.  Tho  advantage 
eight  be  that  the  f lo*  could  look  much  alaplcr  In  Lagranglan  coordinate*. 
The  disadvantage  l*  that  an  Initially  *ell  ordered  ae*h  will  become 
badly  distorted.  The  problem  then  l»  to  keep  track  of  which  partlclea 
are  in  the  immediate  vicinity  of  a  given  particle  In  order  to  determine 
the  force*  they  exert.  Mhtle  work  ha*  oeen  (and  1*  being)  done  on  thl* 
formulation.  It  la  fell  that  the  mean*  to  implement  such  a  calculation 
efficiently  are  not  sufficiently  well  understood  lo  be  used  In  evaluating 
feasibility  of  turbulence  calculation  ho. e  (l.e.,  the  study  of  such 
methods  is  a  worthy  reaearch  project  In  It*  own  right).  Accordingly, 
we  realrict  ouraelve*  here  lo  the  N-S  equation*  in  the  form  given  above. 

(A  recently  proposed  method  of  thorln8*  which  Is  somewhat  related  may 
be  noted.  In  two  dlmenalona  the  method  follow*  the  motion  of  s  finite 
number  of  vortex  element*;  however,  the  accuracy  obtainable  and  the 
extension  lo  three  dimension*  is  not  yet  known.) 

Another  way  of  characterising  computational  methods  Is  according  lo 
whether  they  are  finite  difference  In  ph) steal  space,  or  spectral  methods. 
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or  a  combination. 


The  difference  Method*  can  be  aubdlvidod  into  whether  they  cal¬ 
culate  the  ’prlaltlve"  variable*  1*  and  v  or  related  onca.  (In  tiro 
dimension*  those  are  *  stream  function  and  a  vortlcily  function.  In 
th roe  dimension*  one  need*  a  vector  potential.)  For  t«o-d Isemloml 
calculation*  the  u»o  of  th^ae  derived  variables  has  been  rather  popular. 
However,  there  Is  some  suspicion  that  sore  accurate  results  are  obtain¬ 
able  using  the  prlaltlve  variables.  The  reason  la  that  In  calculating 
the  vortlcity  function  one  Is  trying  to  find  the  derivative  of  a  rapidly 
varying  quantity.  Accordingly  »e  will  make  estimates  assuming  P  and  v 
are  being  calculated  directly. 

In  assessing  finite  difference  net hod*  in  turbulence  calculations 
It  Is  important  to  note  that  turbulence  Is  charact erl red  by  a  wide 
range  of  excited  scales  of  motion,  (See  Section  111  for  sn  example,) 

To  deal  with  this,  one  can  introduce  very  many  grid  points.  Another 
approach,  when  one  knows  a  priori  where  the  small  scale  motions  will 
be,  Is  to  appropriately  stretch  ttc  coordinate  scale.  For  example 
In  the  pipe  flow  (similar  to  the  flow  described  In  Section  11),  most 
of  the  turbulent  production  occurs  at  a  distance  from  the  walls  that 
Is  only  a  ttrv  percent  of  the  radius.  Accordingly,  one  would  like  to 
arrange  the  mesh  so  that  a  significant  fraction  (perhaps  one-fourth) 
is  devoted  to  thut  small  region.  Actually  one  will  probably  want  to 
go  In  both  directions. 

Another  procedure*  is  to  go  to  higher  order  difference  schemes. 

Thus  one  can  apparently  obtain  the  **me  accuracy  (In  d  dimensions)  by 
using  (£)  of  the  mesh  points  needed  In  a  second-order  difference  calcu¬ 
lation  by  going  to  a  fourth-order  scheme.  Tills  can  be  very  significant. 

In  a  two-  or  three-dimensional  calculation,  one  has  to  store  st  sny 
2  3 

lime  step  2 K  or  iS*  numbers,  respectively.  Mncc  the  h‘*s  of  Interest 
are  on  the  order  of  1211  to  102*1,  a  factor  of  four  (or  eight)  Is  of 
considerable  Importance  in  terms  of  memory  requirements.  One,  of  course. 
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pay*  *  price— *ort  arl thmetlcal  operations  **e  needed  and  henco  com¬ 
puting  time  is  increased.  Estimate*  of  shat  on#  gain*  and  loses  by 
going  In  this  direction  are  very  similar  io  the  comparison  given  below 
between  pure  spectral  methods  s  d  a  mined  spoctral-dlf fcrenc#  approach. 

Spectral  methods^  Involve  expanding  the  fields  in  term*  of  an 
appropriate  finite  set  of  ortho  ;ona.i  functions.  Tho  h-S  equations 
then  become  coupled  ordinary  differential  equations  (in  time)  for  the 
expansion  coefficients.  This  formulation  has  a  number  of  attractive 
features.  Thus,  for  a  given  of  accuracy,  one  needs  to  consider 

many  fewer  degrees  of  freedom.  Also  convergence,  *hcn  it  occurs,  Is 
very  rapid  (l.o.,  to  pass  from  ray  55  accuracy  to  15  seems  to  require 
very  few  additional  orthogonal  functions).  The  disadvantages  are  that 
one  must  ;.lnd  functions  with  appropriate  properties  corresponding  to 
a  given  geome  trical  arrangement.  Further  the  transforms  that  the 
expansion  implies  must  be  evaluatable  by  some  efficient  algorithm, 
such  as  the  fast  Fourier  transform.  However,  such  functions  are 
known  for  tho  simple  problems  envisaged  here.  For  periodic  conditions 
one  can  use  Fourier  transforms.  For  pipes  and  channels,  one  can  use 
Chebyahev  polynomials.  Another  disadvantage  Is  an  Increase  In  computing 
time.  (tfc  will  come  back  to  this.) 

Actually,  one  is  not  constrained  to  a  pure  form  of  finite  dif¬ 
ference  methods  or  to  pure  spectral  methods.  A  mix  in  which  some 
coordinates  are  treated  one  way  and  others  another  Is  possible,  snd 
probably  even  optimal. 

To  get  some  possible  estimates  of  computer  requirements,  we 
consider  two  methods  of  calculation.  #*  Probably  neither  is  optimal 
(In  particular,  cf  .  Section  V),  but  they  do  describe  the  kind  of  ranges 
one  might  go  over  if  one  tried  various  of  the  sltcrnatlves  sketched 
above . 
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A .  Method  1 

One  solves  the  velocity  equation  as  a  marching  problem  with  flnlto 

differences  In  apace  and  time  (assuming  a  second-order  scheme).  Assuring 

X  aesh  points  in  each  spatial  direction,  in  two  (three)  dlaenslons  thoro 
2  3 

are  3N  (4N  )  variables  to  be  obtained  at  each  tlae  point.  If  N  Is  large 
the  aaln  computation  problem  la  to  dotoralno  tho  pressure  from  the 
Poisson  equation 

V2 P  -  -c  •  (v  •  t?)v  .  (3) 

Let  us  suppose  that  the  solution  can  be  obtained  by  a  method  equivalent 

to  a  Fourier  transform  technique.  (With  periodic  boundary  conditions, 

this  would  Indeed  bo  a  Fourlor  transform.  For  pipes  and  channels  this 

might  be  a  Fourier  transform  in  two  dimensions  and  Chebyshev  expansion 

In  the  third.)  kc  remark  that  this  is  not  necessarily  the  most  efficient 

method  of  solution,  kith,  for  example,  transforms  In  all  but  one 

dimension  and  treating  tho  reaialnlng  dimension  by  methods  of  inverting 

trldlagonal  matrices  or  Hocknoy's  method''  ono  can  get  by  with  slightly 

fewer  operations.  (Tho  saving  though  Is  only  a  factor  on  the  order  of 

two  or  so,  and  hence  we  lgnoro  It.)  For  large  X,  thon  our  major  computing 

offort  is  to  take  the  Fourlor  transform  of  tho  pressuro  equation,  and 

2  3 

thon  Invert*  This  Involves  tiro  Fourlor  transforms  on  h*  (N  )  variables* 

r  2  2 

Using  the  fast  Fourlor  transform  Algorithm  «ro  then  need  2{N  ln^N 

2  2  3  3 

additions  plus  fN  In^N  multiplications]  or  2fN  In^N  additions  plus 

•  3  3 

JN  lu^N  }•  For  simplicity  lot  us  assume  ore  multiplication  to  be 
equivalent  to  t**>  additions*  (Tho  timer  for  those  is  Indeed  roughly 
that  of  the  Ill lac*) 

Then  the  additions  needed  In  2D  are 

2  2 

4N  In^N  ,  <*U) 
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and  in  3D  wo  need 


4N 


3 


lr. 


X 


3 


(4b) 


for  one  time  step. 

Actually  these  are  only  asymptotic  formulas.  For  the  range  of  N 

wc  will  be  considering  they  are  not  too  bad  for  two-dimensional  estimates, 

but  rather  poor  for  the  three-dimensional  case.  Indeed  we  have  completely 

omitted  the  calculation  of  the  velocity  terras  in  real  space  and  the  Fourier 

transform  of  the  pressure  in  transform  space.  Since  these  are  both  local, 

in  that  at  most  a  point  and  its  close  neighbors  are  involved,  this  givos 

2  3 

a  number  of  operations  ~  X  (X  ).  A  more  correct  estimate  of  the  numbor 
of  equivalent  additions  is  then  in  2D: 

(5a) 

(5b) 


X  2  2  2 

ir  4N  In  X  +  SON 
%/f  2 


per  time  step 


and  in  3D: 


X  3  3  3 

jf  *  4X  In  X  +  150N 


per  time  step 


How  ii*any  time  steps  arc  needed?  For  stability  we  need 


At 


XU 


where  U  is  typical  of  the  overall  flow  velocity  and  L  is  a  typical 
dimension  of  the  system.  The  total  time  T  that  we  might  want  to  cal¬ 
culate  is  some  significant  multiple  (perhaps  about  ten)  of  the  macro- 
scoplc  time  L/U  for  something  to  happen  to  the  whole  system.  The  number 


Note:  The  fact  that  we  may  need  a  fine  grid  resolution  near  a  wall  may 
still  not  cause  trouble  with  the  Courant  stability  condition.  Thus  in 
the  boundary  layer  the  flow  is  nearly  parallel  to  the  wall  and  the  condi¬ 
tion  involves  only  the  projection  of  tho  convecting  velocity  on  the  grid 
increment  in  that  direction. 
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of  time  steps  N  then  is  such  that 
T 

N  At  =*  mL/U  or  N  **  mN 

T  T 

i 

i 

(6a) 

(6b) 

The  total  time  for  such  a  calculation  is  then  the  above  numbers  times 
t,  the  addition  time. 

B.  Method  2 

Here  one  imagines  using  purely  spectral  methods.  For  example, 
suppose  we  can  use  Fourier  transforms.  If  k  is  the  maximum  wave  vector 
in  each  direction  we  have 

2 

in  two  dimensions  2  *  (2k)  real  variables 

3 

and  in  three  dimensions  3  *  (2k)  real  variables  . 

Here  the  whole  problem  is  involved  with  calculation  of  the  nonlinear 
terms  in  the  velocity  equation  (which  were  asymptotically  unimportant  in 
Method  1).  An  efficient  means  of  doing  this  is  by  passing  back  and  forth 
between  Fourier  space  and  real  space.  It  has  been  shown  that  in  two 
(three)  dimensions  this  can  be  done  with  20  (72)  fast  Fourier  transforms. 

Therefore,  for  one  time  step  the  number  of  effective  additions  is 

k  2  2 

Jf  ~  40<2k)  In  (2k)  in  two  dimensions 

6/7  a  2 

and 

k  3  3 

Jf  ***  144(2k)  In  (2k)  in  three  dimensions, 

a  » 
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Therefore  the  total  number  of  operations  is 


3  9  3  2 

in  2D  =  mN  {41n  N~  +  50}  **  N  {401nrN  +  500} 

2 


4  3  4  3 

and  in  3D  =  mN  (41n  N  *  +  150)  »  N  {401n  N  +  1500) 

2  2 


-  -  _ ■ _ _ _ _ 


of  time  steps  Nt  then  is  such  that 

N  At  =“  mL/U  or  N  »  mN 

i  T 

Therefore  the  total  number  of  operations  is 

i n  20  =  mN3{41n2N2  +  50}  »  N3{401n2N2  +  500} 

and  in  3D  =  mN4{41n2N3  +  150}  N4{401n2N3  +  1500} 

The  total  time  for  such  a  calculation  is  then  the  above  numbers  times 
t,  the  addition  time, 

B.  Method  2 

Here  one  imagines  using  purely  spectral  methods.  For  example, 
suppose  we  can  use  Fourier  transforms.  If  k  is  the  maximum  wave  vector 
in  each  direction  we  have 

in  two  dimensions  2  •  (2k)2  real  variables 

and  in  three  dimensions  3  •  (2k)3  real  variables  . 

Here  the  whole  problem  is  involved  with  calculation  of  the  nonlinear 
terms  in  the  velocity  equation  (which  were  asymptotically  unimportant  in 
Method  1).  An  efficient  means  of  doing  this  is  by  passing  back  and  forth 
between  Fourier  space  and  real  space.  It  has  been  shown27  that  in  two 
(three)  dimensions  this  can  be  done  with  20  (72)  fast  Fourier  transforms. 
Therefore,  for  one  time  step  the  number  of  effective  additions  is 

k  2  2 

Jf  ~  40(2k)  In  (2k)  in  two  dimensions 

and 

144 (2k)  ln^(2k)  in  three  dimensions. 


(6a) 

(6b) 
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For  comparison,  it  is  important  to  note  that  the  accuracy  for  a 
given  k  is  that  by  second  order  finite  differences  for  an  equivalent 
mesh  of  size 

N  =  4k  . 
e 

Expressing  everything  in  terms  of  we  see  that  the  spectral  method 
requires 


variables  in  two  dimensions 

variables  in  three  dimensions, 

2  3 

This  can  be  compared  with  the  numbers  2N  and  4N  for  Method  1  (i.e,, 
in  two  dimensions  we  have  a  saving  in  memory  by  a  factor  of  four  and 
in  three  dimensions  by  more  than  a  factor  of  10). 

What  do  we  pay  in  computing  time?  Using  the  asymptotic  formulas 
we  obtain 

in  two  dimensions 

and  in  three  dimensions 

Actually  the  time  disadvantage  here  is  quite  overestimated.  The 
relatively  s...all  N  that  would  be  of  interest  is  such  that  these  ratios 
are  more  like  a  factor  between  one  and  two.  The  main  disadvantage  of 
the  spectral  methods  seems  more  likely  to  be  their  lack  of  generality 
and  increased  complexity  of  coding. 


Jf 

Jr 


Jr 

jc 


=  2.5 
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a  3 
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V  ILLIAC  IV  CHARACTERISTICS 


To  estimate  computing  times  it  is  useful  to  have  a  brief  description 
of  the  properties  of  Illiac  IV,  as  perhaps  it  may  be  within  a  year. 

(Only  characteristics  essential  for  our  evaluation  are  given.)  The 
essential  conclusion  is  that  memory  requirements  drastically  limit 
the  usefulness  of  the  Illiac  in  turbulence  calculations.  The  essential 
characteristics  include : 

(1)  Processors .  There  are  64  parallel  processors.  These  can 
all  be  performing  the  same  operations,  or  some  can  be 
inoperative.  One  PE  can  do  a  64-bit  floating-point  add 
in  30  ns  and  a  floating-point  multiply  in  600  ns. 

(2)  Memory . 

(a^  There  is  a  1.28  X  10  64 -bit -word  working-array  storage. 

g 

(b)  There  is  a  16  X  10  64-bit-word  disk  storage.  The 
transfer  rate  from  disk  to  array  is  10^-word  per  second. 
This  is  somewhat  misleading,  so  we  spell  out  the  details. 
The  disk  is  arranged  in  52  bands,  each  with  300  pages 

of  10^  64-bit  words.  The  unit  of  transfer  is  one  page. 
Then  133  /is  are  needed  to  transfer  a  page.  However, 
this  requires  having  immediate  access  to  the  desired 
page.  The  disk  rotates  in  40  ms.  If  one  has  to  seek 
out  a  random  page  we  might  expect  the  time  to  be  -10  to 
20  ms.  Further  limitations  are:  (i)  it  takes  133  us 
to  switch  between  bands,  and  (ii)  one  cannot  read  on 
one  band  and  write  on  another  simultaneously .  Indeed 
it  takes  another  133  jus  to  change  from  reading  to 
rewrite.  The  principal  conclusion  here  is  that  to 
use  the  10  -words-per-second  transfer  rate  one  must 
have  stored  data  in  exactly  the  right  place.  As  much 
as  Dossihle  all  numbers  needed  during  a  set  of  calcula¬ 
tions  should  have  come  from  a  given  page  and  neighboring 
pages . 

(c)  There  is  also  a  buffer  disk  of  abcut  ten  times  the  Illiac 
disk  and  a  transfer  rate  about  one-tenth  that  of  the 
Illiac  disk. 
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(d)  There  is  a  Unicon  laser  storage  device  with  an  on-line 
storage  capacity  of  ~109  64-bit  words  and  infinite  off¬ 
line  storage  with  a  transfer  rate  of  105-words-per- 
second . 

Properties  (c)  and  (d)  have  little  effect  on  computing  capability 
but  do  restrict  the  amount  of  information  one  might  want  to  bring  out 
of  the  computer. 

It  is  somewhat  instructive  to  compare  the  I  Iliac  IV  with  other 
machines.  Take  as  a  basis  the  roughly  comparable  IBM  360/195  and  the 
CDC  7600.  Their  characteristics  are  approximately: 

•  Add  time  ~  100  ns 

•  Multiply  time  ~  200  ns 

g 

•  Fast  Memory  (0.5  to  1)  X  10  14-bit  words. 

We  note  that  since  the  arithmetic  times  are  of  the  order  one-third 

that  of  the  Illiac,  the  factor  of  64  (due  to  the  many  processors)  is 

effectively  reduced  to  a  factor  20.  There  is,  however,  the  question  of 

the  efficiency  with  which  one  can  employ  the  64-fold  parallelism.  For 

3  6 

example,  suppose  one  were  calculating  using  a  (128)  ~  10  grid.  Let 

us  suppose  that  the  bulk  of  the  calculation  can  be  done  with  the  full 

2 

64-fold  parallelism.  The  boundary  points  are  of  order  6  X  (128) 

~  60,000  points.  These  are  about  10%  of  the  calculation  and  might 
require  only  six  processors.  This  will  then  double  the  computing  time. 

A  not  unreasonable  estimate  of  the  Illiac  advantage  over  the  IBM  360/195 
or  CDC  7600  is  then  about  ten.  (This  assumes  the  whole  problem  be 
contained  in  the  array  memory:  If  not,  for  the  reasons  given  above, 
this  advantage  may  be  drastically  reduced.) 

Similar  estimates  for  Star  and  ASC  are  more  difficult  since  much 
less  is  known  about  how  well  these  will  work.  Guesses  as  to  the 
improvement  of  Star  over  the  CDC  7600  seem  to  be  a  factor  of  three 
while  the  ASC  improvement  might  approximate  two  (for  one  pipe)  to 
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eight  for  four  pipes.  Again  the  problem  of  implementing  the  parallelism 


remains  a  question.  (One  interesting  property  of  the  ASC  is  the  apparent 

0 

capability  of  increasing  the  fast  memory  to  (4  to  16)  X  10  words.  This 
could  be  a  great  advantage  for  turbulence  calculations.  ) 

Let  us  now  try  to  make  some  estimates  of  computing  times  for  problems 
of  various  sizes. 

An  absolute  minimum  is  obtained  by  taking  the  number  of  equivalent 

additions,  multiplying  by  the  add  time  and  assuming  we  can  use  the  full 

64-parallelism  perfectly.  In  Tables  2  and  3  we  give  the  results  in  two 

and  three  dimensions  corresponding  to  Methods  1  and  2.  We  give  number 

of  variables,  time  per  step,  and  total  time  (assuming  T  =10NAt). 

(Here  t  ,  t  are  time  per  step  for  Methods  1  and  2;  U  ,  U  the  number 
y  c  12 

of  variables;  and  T  ,  T  are  the  corresponding  total  times.)  An  impor- 

-L  6 

tant  point  to  be  noted  from  these  tables  is  that  for  comparable  accuracy 
t2  is  not  much  greater  than  t  ,  while  the  number  of  variables  to  be  ke,:t 
is  much  smaller. 

To  be  realistic  we  should  allow  some  factor  for  inefficiency. 
Accordingly,  all  times  should  really  be  multiplied  by  some  factor-- 
probably  between  two  and  ten:  Let  us  say  five.  What  problems  are 
now  within  reasonable  limits?  Suppose  we  put  an  upper  limit  of  100 
hours  for  one  problem.  Then  we  see  that  in  two  dimensions  we  can  at 
most  consider  N  =  2048  (or  K  «  512),  while  in  three  dimensions  the 
limits  are  N  =  256  (K  =  64).  If  we  drop  our  limit  to  approximately 
ten  hours  the  numbers  are  N  =  1024  (K  =  256)  in  two  dimensions  and 
N  =  128  (K  =  32)  in  three  dimensions.  We  note  that  while  these 
estimates  are  quite  rough,  it  is  clear  that  the  rapid  rise  in  total 
computing  time  with  increasing  N  suggests  these  upper  limits  cannot 
be  wrong  by  more  than  a  factor  of  two. 
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However,  in  considering  using  i'll  Inc  IV  for  such  turbulence  cal¬ 
culations  on©  should  also  consider  the  relatively  small  array  memory* 

We  describe  the  three-dimensional  case  in  the  following  paragraph. 

Suppose  one  restricts  attention  to  problems  that  can  be  done  using 
only  the  fast  array  memory.  Presumably  at  any  time  wo  need  to  be  storing 
-2V  words .  For  Method  1  wo  aro  restricted  to  X  <  32  and  for  Method  2 
to  K  <  1G.  Clearly,  In  order  to  soivo  interesting  problems,  use  will 
have  to  be  made  of  the  I  lilac  disk.  The  problems  that  are  then  possible 
are  thoae  with  accuracy  comparable  to  that  achievable  with  second-order 
*  ini  to  differences  with  128  -SC  256.  Methods  like  Type  1  are  clearly 
limited  by  disk  space  to  X  128.  Here  the  disk  is  completely  filled 
(Indeed  slightly  over)  by  numbers  one  needs  to  be  keeping.  The  situation 
Is  similar  in  methods  of  Type  2  for  K  •  64.  In  any  case  one  Is  faced 
with  a  situation  where,  at  any  lnstanco,  almost  all  the  numbers  are 
being  stored  in  the  disk.  Then  the  disk  transfer  rate  also  gives 
some  time  llmitstlons. 

Thus,  consider  a  situation  such  as  Method  1  with  X  128.  (We 

say  approximately  alnce  the  exact  caae  does  not  fit  in;  however,  w© 

Imagine  this  taken  care  of,  for  example,  by  going  to  32-blt  words.) 

At  any  time  stop  we  are  storing  -6  x  10°  words  from  a  previous  time 
_  0 

step.  The  remaining  8  X  10  memory*  places  are  to  be  repeatedly  used 

during  the  time  step.  If  coding  could  br  done  perfectly,  one  could 

Imagine  proceeding  so;  The  6  x  10  velocity  variables  from  a  given 

2 

step  arc  fed  in.  From  those  the  torn*  v  •  tTv,  V  v  and  V  •  (v  ■  ?v) 
are  computed.  Since  these  are  so  numerous  they  must  be  taken  out  to 
the  disk.  Fourier  transforms  are  then  taken  plane  by  plane— again 
they  must  be  taken  out.  Next  one  transforms  In  the  direction  perpen¬ 
dicular  to  the  plancs  —  thls  la  again  a  pass  in  ami  out.  The  same  must 
be  done  to  invert  the  transform  of  the  solution  of  the  Poisson  equation. 
Finally,  all  torms  must  be  combined  to  obtain  the  next  values  of  the 
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velocities.  Thu*  there  sr*  *  number  of  passe*  (about  ten)  In  and  out 
of  the  disk,  (.Sot  all  the  passes  require  a  full  8  x  10**  words.)  How¬ 
ever,  the  not  effect  is  that  something  like  ft  y  (8  x  106)  wo.'ds  have 
to  go  in  and  out  of  the  disk.  With  a  transfer  rate  of  10^  words  per 
second,  this  is  on  the  order  of  six  second*.  We  remark  thai  with  our 
previous  eat  last*  of  (5  X  2.3)  seconds  for  computing  time  pe.-  step, 
this  does  not  Materially  change  the  tiao  that  the  total  calculation 
would  take.  It  does  indicate  that  just  enlarging  tne  disk  wan  fold 
does  not  change  tho  situation  as  to  what  probice*  can  *»#  attacked 

The  real  lialtatlon  lapoacd  here  by  tho  disk  Is  hidden  tv  our 
assuapt Ion  of  perfect  coding.  If  at  any  given  pass  all  the  numbers 
are  not  on  the  a.ute  page  with  all  those  needed  to  coapute  th«  next 
stage  fro*  thesi,  and  Indeed  all  pages  are  not  arranged  to  be  accessible 
In  sequential  order,  our  assumed  transfer  of  rate  of  103  words  In 
133  us  drops  to  10  as  or  less.  Clearly,  to  perfora  calculations  el  ere 
s  large  fraction  of  the  disk  mist  be  repeatedly  read  and  rowrlt'cn, 
we  require  s  very  careful  choice  of  method  of  calculation  and  greet 
care  In  coding.  It  would  seea  that  the  aaln  determinant  for  the 
method  of  calculation  may  well  be  the  eaae  with  which  It  permits  coding 
for  efficient  access  and  egroaa  from  the  dlak.  Thla  la  a  difficult 
problea  and  will  take  careful  atudy.  However,  we  bellove  It  not 
necesaarlly  Insoluble.  (In  Appendix  0,  a  possible,  though  probably 
not  opt  las  1 ,  approach  la  sketched.) 
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velocities.  Thus  there  are  a  number  of  passes  (about  ten)  In  and  out 
of  the  disk.  (f*'ot  all  the  passes  require  a  full  8  x  10  words.)  How¬ 
ever.  the  not  effoct  Is  that  something  llko  5  X  (8  X  10®)  words  hsvc 
to  go  in  and  out  of  the  disk.  With  a  transfor  rate  of  10?  words  per 
socond,  this  is  on  the  order  of  six  seconds.  We  remark  that  with  our 
previous  estimate  of  (3  X  2.3)  ceconds  for  computing  time  per  stop, 
this  does  not  materially  chango  the  time  that  the  total  calculation 
would  take.  It  does  Indicate  that  Just  enlarging  the  disk  manyfold 
does  not  change  the  situation  as  to  what  problems  can  be  attacked. 

The  real  limitation  Imposed  here  by  the  disk  Is  hidden  by  our 
assumption  of  perfect  coding.  If  at  any  given  pass  all  the  numbers 
are  not  on  the  aamo  page  with  all  those  needed  to  compute  the  next 
stage  from  them,  and  lndoed  all  pages  are  not  arranged  fo  be  accessible 
In  sequential  order,  our  assumed  transfor  of  rate  of  103  words  in 
133  us  drops  to  10  as  or  leas.  Clearly,  to  perform  calculations  where 
a  large  frnctlon  of  the  disk  must  be  repeatedly  read  and  rewritten, 
we  require  a  very  careful  choice  of  method  of  calculation  and  groat 
c-'re  In  coding.  It  would  seem  that  the  main  determinant  for  the 
method  of  calculation  may  well  be  the  case  with  which  It  permits  coding 
for  efficient  access  and  egross  from  the  disk.  This  Is  a  difficult 
problem  and  will  tako  careful  study.  However,  we  believe  It  not 
necossarlly  Insoluble.  (In  Appendix  D,  a  possible,  though  probably 
not  optimal,  approach  Is  sketched.) 
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VI  CONCLUSIONS 


In  principle  the  advent  of  a  new  generation  of  computers  (in  parti¬ 
cular  the  Illiac  IV)  Makes  possible  the  numerical  solution  of  some 
interesting  problems  in  turbulence.  By  themselves,  however,  they  will 
not  give  the  ideal  goal,  an  ability  to  predict  the  flow  for  a  wide 
variety  of  practical  problems.  Uany  basic  questions  and  problems 
involve  Reynolds  numbers  far  beyond  present  capabilities.  What  can 
be  achieved  are: 

(1)  The  obtaining  of  solutions  to  a  number  of  benchmark  problems. 
These  can  be  used  to  test  approximation  methods  for  Reynolds 
numbers  R  <  1C"*  -  4  x  10^ . 

(2)  The  obtaining  of  greater  insight  into  the  nature  of  turbulence. 
This  would  be  the  result  of  a  number  of  numerical  experiments. 

Such  an  effort  seems  desirable;  probably  i'.  does  not  require  a 
massive  program. 

Doing  these  problems  on  the  Illiac  IV  presents  severe  difficulties. 
The  army  memory  is  small  compared  to  the  number  of  independent  variables 
needed,  at  least  for  some  of  the  more  obvious  methods  of  computation. 

If  desired,  a  reasonable  program  might  be  the  following: 

(1)  An  attempt  should  be  made  to  see  how  a  typical  problem  could 

be  coded  so  ns  to  achieve  a  resolution  of  distances  of  the 

order  of  1/128  of  the  overall  scale  (e.g.,  in  a  pipe  one  wants 

to  resolve  distances  -  1/128  of  r  ).  The  method  of  calculation 

o 

used  will  bo  primarily  determined  by  the  size  of  the  array 
memory.  One  approach  is  to  see  how  one  can  program  so  as  to 
efficiently  access  the  disk  memory.  If  it  is  found  that  this 
cannot  be  done  well,  an  alternative  might  be  to  go  to  very 
high-order  difference  schemes.  Presumably  a  significant 
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reduction  of  the  number  of  independent  variables  can  be 
achieved,  at  the  expense  of  an  increase  in  computing  time, 
Tho  exnct  advantngos  and  disadvantages  obtained  in  this 
wny  are  not  known  and  would  have  to  be  dotormined. 

(2)  Assuming  this  coding  problem  has  been  solved,  one  would 
calculate  tho  problems  sketched  in  Section  III.  Sinco 
many  of  the  results  are  desired  for  theoretical  insight, 
the  calculation  should  be  carefully  followed  by  on  advisory 
i’roup.  The  actual  choice  of  what  to  calculate  next  will 
Jepend  on  what  has  alroady  boon  found. 

Since  huge  quantities  of  numbers  will  bo  generated, 
it  will  be  nocossary  to  be  rather  selective  as  to  what  one 
take-  out  of  the  computer.  For  example,  one  might  want 
only  to  read  out  every  10  to  20  time  steps  the  various 
quantities  for  which  experimental  results  are  available 
(see  Section  II-B).  In  particular,  though,  one  would  want 
the  energy  spectrum,  vorticity  spectrum,  skewnoss,  and 
kurtosis . 

(3)  Since  the  amount  of  experimental  Information  at  present 
far  exceeds  the  amount  of  numerical  results,  an  experi¬ 
mental  program  does  not  now  seem  in  order.  When  (if?) 
numerical  results  are  available,  there  may  be  some 
questionable  features  which  would  suggest  specific 
experiments. 

(4)  Even  if  tho  problem  of  appropriate  coding  for  the  llliac  IV 
proves  Insoluble,  some  parts  of  tho  above  program  could  be 
usefully  done  on  a  lessor  computer.  For  example,  even  wit), 
nn  IBM  360/1 ^5  or  a  CDC  7600  some  of  tho  above  problems 
could  be  done,  but  at  Reynolds  numbers  lower  than  those 
envisaged  for  the  llliac. 

It  is,  of  course,  possible  thnt  the  computing  estimates  thnt  have 
been  made  are  too  high,  in  thnt  alternative  methods  of  calculation  mny 
require  less  computing  time.  However,  It  Is  thought  thnt  our  estimates 
nre  fairly  conservative.  In  any  event  the  primary  conclusion — thnt 
the  computations  needed  Increase  dramatically  with  Reynold’s  number— 
is  certainly  correct. 
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A  NAVIER-STOKES  COMPUTER? 

1 .  Summary 

The  conceptual  design  of  a  special-purpose  computer  array  to  solve 

the  N-S  (Navier-Stokes)  equation  is  discussed  below.  An  array  of  100 

to  150  special  processors  (each  roughly  like  a  bare  modern  minicomputer) 

with  a  100-ns  add  time  will  permit  the  solution  of  the  second-ordei 

difference  approximation  to  the  N-S  equati  n  in  a  total  processing  time 
4 

(for  10  time  steps)  of  several  minutes  (N  =  64)  to  seven  hours  (N  =  1C24), 

This  processing  time  could  be  decreased  by  paralleling  and/or  pipelining 

5 

in  each  processor.  The  total  hardware  cost  is  of  the  order  of  $5  X  10 
6 

to  10  . 

The  major  problem  is  not  the  processors  but  memory.  A  total  of 
7  11 

6  X  10  (N  =  64)  to  2  X  10  (N  =  1024)  bits  are  required.  It  appears 

-4 

that  for  the  largest  calculation,  a  cost  of  no  more  than  $10  per  bit 

is  required.  In  order  to  match  memory  speed  to  processor  speed,  each 

processor  will  require  16  K  to  32  K  of  100-ns  memory  in  at  least  eight 

-3 

banks  and  four  or  eight  drums  or  disks  with  about  an  8  X  10  -s  average 

—6  6 
access  time,  a  5  X  10  -s-per-word  transfer  time  and  1  to  2  X  10  total 

6 

storage.  That  is,  a  total  of  1.6  to  3.2  X  10  words  of  100-ns  memory 

6 

and  400  to  800  independent  disks,  each  with  a  capacity  of  1  to  2  X  10 

words  is  required  to  handle  even  the  N  =  512  problem.  The  total  cost 

7 

of  this  memory  will  certainly  be  in  excess  of  $10  , 
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2. 


Discussion 


Assume  that  there  are  K  V-Boxes,  L  F -Boxes,  and  M  P-Boxes,  and  a 
memory  (all  are  described  in  the  following  sections),  then  the  total 
time  to  calculate  V  and  P  on  all  mesh  points  for  one  time  step  is 


/T  \  /T  \  /T  \ 
(  _V  1  |  _F  j  [  JP  \ 

Vk/Vl/Vm/ 


=  N 


3 


The  total  cost,  C^,  is  the  sum  of  the  cost  of  the  V,  F,  and  P  boxes, 
the  C  box,  and  memory,  and  is 


C  =  KC  +  LC  +  MC  +  C  +  C 
T  V  F  P  c  M 


=  (K  +  L  +  M)(c  +  c  +  8c  +  c)+C  +C 
am  r  w  c  M 


where 


c  =  Cost  of  one  adder 
a 


c  =  Cost  of  one  multiplier 
m 


c  =  Cost  of  one  register 


c  =  Cost  of  wiring, 
w 

The  memory  cost,  C  ,  is 
M 


C  =  224  c  N3,  or 
M  o 
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°M  =  64  V  +  32  c1(7n3  “  2n) 

=  64  (c  -  c  )n  +  224  c  N3 
o  1  1 

depending  on  whether  a  single  high-speed  memory,  with  a  cost  per  bit 
of  Cq  or  a  two-level  memory  (2n  words  of  fast  memory  at  a  cost  per  bit 
of  Cq  and  the  remainder  slow  memory  at  a  cost  per  bit  of  Cj)  is  used. 
The  high-speed  memory  has  a  cycle  time  t  St  and  the  slow-speed 
memory  has  a  cycle  time  of  t  . 

Assume,  for  the  sake  of  making  an  estimate,  that  roughly  the  same 
tj<ne  is  spent  in  each  of  the  boxes.  That  is, 


-  ii  R3 -  (=55  - 

K  L  M 

Taking  N  =  256,  it  is  seen  that  M  =  K/10  and  L  =  K/4  will  satisfy. 

Then, 

T  «  (600  N3/K)  t 

O 

Taking  N  in  the  range  64  to  1024  and  tQ  =  100  ns  (currently  available 
adders  achieve  this),  it  is  seen  that  K  -  100  gives  T,  the  total  compute 
time  per  time  step,  in  the  range  lo"2  to  60  s,  for  a  maximum  of,  say, 
seven  hours  per  computation. 

If  seven  hours  per  computation  is  too  long,  it  can  be  considerably 
shortened,  by  at  least  an  order  of  magnitude.  If  the  (second  order) 
difference  approximations  to  the  Navier-Stokes  equations  are  examir  d 
in  detail,  it  becomes  clear  that  it  is  possible  to  build  into  each  of 
the  processors  a  fairly  high  degree  of  parallelism  and/or  pipelining 
and  thus  achieve  the  speedup.  However,  there  seems  to  be  little  point 
in  this  because,  as  will  be  seen,  the  memory  requirements  are  already 
difficult  and  speeding  up  the  processors  will  only  make  them  worse. 
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I 

In  order  to  keep  the  processors  busy  it  is  necessary  to  supply 

some  memory  with  an  access  time  t  , 

o 

t  £  t  =  100  ns, 
o  o 

If  w*)  had  a  set  of  main  memories,  one  for  each  processor,  it  is  clear 

0 

that  even  for  an  N  =  64  calculation  something  of  the  order  of  2  X  10 
6 

words  (64  X  10  bits)  of  memory  would  be  required.  If  100-ns  memory 

5 

costs  $0.01  per  bit,  the  cost  is  of  the  order  of  $6  X  10  .  For  N  =  1024, 

9 

the  co5> t  is  of  the  order  of  $2  X  10  . 

If  a  number  of  secondary,  slow  speed,  memories  are  used,  then  the 

access  time  (see  below)  could  be  in  the  range  1  to  10  jus  depending  on 

the  number  of  banks  (one  to  eight)  in  each  of  secondary  memories.  Still, 

7  11 

between  6  x  10  and  2  X  10  bits  are  required.  If  the  cost  per  bit 

-4  0 

were  only  $10  the  largest  calculations  require  of  the  order  of  $20  X  10 
worth  of  memory. 

Going  to  drums  or  disks  for  the  secondary  memory  will  help,  but  not 

too  much.  If  each  primary  memory  (2n  words,  100-<us  access  time)  has 

l  disks,  each  with  an  average  access  time  of  t  and  a  transfer  time  per 

o 

word  of  t  ,  then 

T  T 

a  l  <  n  -6 

—  +  —  ^  2  X  10  s 
n  l 

The  right-hand  side  is  an  estimate  of  the  time  the  fastest  box  (V,  F, 
or  P)  takes  to  penorm  its  calculation  on  one  mesh  point.  Taking 

-3 

t  =  8  X  10  s 
o 

-6 

T  =  5  x  10  s 
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A  number  of  pairs  of  n  and  &  can  be  chosen  to  satisfy  this  inequality. 

Some  possible  choices  are 

n  =  8  K,  Si  =  4 

n  =  8  K,  H  =  8 

n  =  16  K,  i  =  4. 

The  first  gives  a  time  just  a  bit  too  long  and  the  others  satisfy  the 
ineaual ity  comfortably.  Therefore,  100  fast  memories  (100  ns)  either 
of  16  K  or  32  K  words  capacity  and  400  or  800  independent  disks  .  each 

with  10  to  2  X  10  words,  are  required  to  handle  the  N  »  1024  problem. 

4  6 
If  each  disk-drive  costs  only  $10  ,  the  disk  cost  is  $4  to  8  X  10  . 

7  8 

Also  5  X  10  to  10  bits  of  fast  (100  ns)  memory  is  required.  If  the 

-3  6 
cost  were  only  $10  per  bit  the  cost  is  $5  to  10  x  10  . 

Finally  a  word  on  the  processor  costs.  Each  of  these  boxes  appears 

roughly  equivalent  to  a  good  minicomputer.  The  bare  cost  of  the  processor 

3 

should  be  similar,  say  $5  to  10  X  10  .  This  is  equivalent  to  taking 

c  ~  -be  ^  c  ^  $100  to  $200 
a  m  r 

and  assuming  that  the  wiring  cost  is  two  to  three  times  the  hardware 

5  6 

cost.  The  total  processor  cost  is  then  of  the  order  $5  X  10  to  10  , 
which  is  considerably  smaller  than  the  memory  cost. 

3 .  The  V  Box 

The  V  box  calculated  V  and  the  right-hand  side  of  the  Poisson 

equation  for  P  on  one  mesh  point.  Assume  that  all  of  the  operations 

necessary  to  perform  these  calculations  are  performed  sequentially  and 

that  the  adder  and  multiplier  are  not  pipelined.  The  total  time  for 

this  calculation,  t  ,  is  then 

v 


t  =  (number  of  addition) (addition  time) 
v 

+  (number  of  multiplications) (multiplication  time). 

These  calculations  require  approximately  100  additions  and  50  multipli¬ 
cations,  Therefore, 


t  =  200  t 
v  o 

3 

The  total  time  to  calculate  on  N  points,  assuming  no  parallel  V  boxes, 
is 

T  =  (200  N3)  t 
V  o 

In  the  wors*  case,  two  successive  operations  will  require  four 
different  operands.  To  „void  waiting  for  the  second  pair  of  operands 
to  be  fetched  from  memory,  they  should  be  fetched  while  the  first 
operation  is  being  performed.  This  ) equires  two  registers  to  hold 
the  current  operands  and  two  registers  *o  hold  the  next  operands.  A 
result  register  to  hold  the  result  of  the  current  operation  is  needed 
as  well  as  a  holding  register  to  save  the  partial  sum  accumulated  in 
calculating  V,  Finally  two  registers  for  constants  are  probably  needed. 
In  total: 

•  Two  current  oporand  registers 

•  Two  next  operand  registers 

•  Two  constant  registers 

•  One  result  register 

•  One  holding  register 

for  a  total  of  eight  registers  are  needed.  In  addition,  one  adder  and 
one  multiplier  are  also  required. 
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The  cost  of  the  hardware  in  a  V  box,  and  assembly  cost,  is 

C  =  c  +  c  +  8c  +  c 
V  a  m  r  w 

where 

Cy  =  Total  cost  of  one  V  box 

c  -  Cost  of  one  adder 
a 

c  =  Cost  of  one  multiplier 

=  Cost  of  one  register 

c  =  Cost  of  wiring  up  a  V  box. 
w 

4 .  The  F  Box 

n 

The  F  box  is  just  a  machine  to  calculate  the  FFT  on  N  points  in 
a  fixed  plane  and  then  over  all  planes,,  A  large  number  of  FFT  machines 
have  been  built  and  it  is  probably  more  accurate  to  estimate  the  time, 
Tf>  to  calculate  the  required  FFT's  as  well  as  the  cost,  C  ,  of  a  single 
FFT  machine,  oy  using  the  speed  and  cost  figures  for  existing  machines, 
than  it  is  to  make  an  a  priori  estimate.  However,  for  completeness, 
such  an  a  priori  estimate  will  be  made. 

In  order  to  calculate  the  FFT  in  one  plane  (NXN)  requires 
2  2 

N  ln^N  additions  and, 

2  2 

ln^N  multiplication,  or 
2  2 

2N  ln^N  operations. 

This  must  be  done  for  N  planes,  requiring 
3  2 

2N  ln^N  operations. 
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After  solving  the  resultant  set  of  tri-diagonal  equations  in  trans¬ 
form  space,  the  inverse  must  be  calculated,  just  doubling  the  number  of 
operations . 

Therefore  the  time  to  take  all  the  FFT  and  inverses,  assuming  no 
parallel  F  boxes,  is 

T  =  (4N3ln  N2)  t 
F  2  o 

It  appears  that  an  F  box  would  have  an  architecture  similar  to  a  V  box, 
w<o,  for  purposes  of  estimate,  it  will  be  assumed  that  their  costs  are 
the  same — i  .e. , 


C 

F 


5 .  The  P  Pox 

The  P  box  solves  the  set  of  tri-diagonal  equations  in  transform 

2 

space.  There  are  a  total  of  N  sets  of  tri-diagonal  equations  to  be 

solved.  It  can  be  shown  that  each  set  of  equations  requires  about  13N 

3 

operations  and  thus  the  complete  solution  requires  13N  operations. 

The  total  time  to  do  this  calculation,  again  assuming  r, :>  parallel  P 
boxes,  is 

T  =  (13N3)  t 
P  o 

Again,  it  seems  reasonable  to  assume  the  cost  of  the  P  box  is 
about  the  same  as  the  cost  of  the  V  box,  so 


C 

P 
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6. 


Tho  C  Box 


The  C  box  is  a  control  unit  for  tho  V,  F,  and  P  boxes*  It  would 

be  a  "hardwired’1  control  unit  which  would  normally  sequence  the  boxes 

to  execute  tho  algorithms  which  calculate  and  P.  It  would  certainly 

be  of  complex  design,  but  probably  no  more  so  than  the  control  unit 

for  any  large  genoral-purpose  computer.  Its  operations  would,  of  course, 

be  concurrent  with  the  other  boxes  and  would  cost  nothing  in  processing 

time.  Le"  C  be  its  construction  cost, 
c 


7 .  Memory 


There  are  a  total  of  N  mesh  points.  A  value  of  and  P  must  be 

3 

stored  for  each  mesh  point;  thus  at  least  4N  words  of  memory  are  re- 

•4 

quired.  If  a  centered-time-dif f erence  scheme  is  used,  then  V  at  two 

3 

time  steps  must  be  stored,  raising  the  total  to  7N  words  of  memory. 

3 

It  is  easily  seen  that  another  3N  words  is  not  required  to  store  the 

new  values  of  V.  Depending  on  the  spatial  differencing  scheme,  only 
2 

a  multiple  of  N  additional  words  is  needed.  For  example,  if  a  second- 

•4 

order  scheme  is  used,  only  the  values  of  V  and  P  in  the  plane  above  and 

■"4 

below  are  needed  to  calculate  V  in  a  particular  plane.  When  the  calcula- 

-*n+l  -»-nn+1  -^_1 

tion  of  V  and  {7  •  L(V  •  7)VJJ  for  plane  2  are  complete,  V  and 

pn  for  plane  1  will  never  again  be  needed.  Therefore  the  values  of 
.+n+l  ^n-1 


V  ,  for  plane  1,  can  be  stored  in  the  memory  holding  V  of  plane  1 
f  r  v  -4  ,  n+1 

and  the  values  of  {7  -  L(V  •  7)VJj  of  plane  1  can  be  stored  in  the 

memory  holding  pn  of  plane  1.  Thus  it  is  necessary  only  to  carry  along 

2  _*H+1  -+  —4  n+l  2 

two  planes  (2N  )  of  V  and  [7  •  [(V  •  7)V]}  for  a  total  of  8N 

additional  words  of  memory.  Because  N  »  8  in  calculations  of  interest, 

3 

this  will  be  neglected  in  comparison  to  7N  . 
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Total  memory  requirements  are  then 

3 

7N  words, 

and,  assuming  each  word  has  32  bits,  which  seems  to  be  about  the 
minimum  acceptable,  a  total  of 

22‘1N3  bits 

is  required. 


To  permit  the  V,  F,  and  F  boxes  to  operate  without  being  held 
up  waiting  for  memory  fetches,  it  is  necessary  to  have  two  words  be 
fetched  from  memory  in  one  add  timo.  Therefore, 


t  *  it 
o  o 


where  is  the  memory  cycle  time, 

3 

It  may  not  be  feasible  (due  to  cost)  to  provide  7N  words  of 

fast  memory.  If  so,  a  buffer  or  cache  memory  of,  say,  2n  words  with 

a  cycle  time  of  t  ,  and  a  secondary  or  main  memory  of  m  words  with  a 
o 

longer  cycle  time  of  can  be  used.  Clearly 


2n  +  m  s= 


Pages  of  n  words  will  be  exchanged  between  the  buffer  and  the 

main  memories.  Because  t  <  t  ,  the  time  required  to  effect  the 

o  1 

exchange  will  be  determined  by  and  will  be  Since  the  pages 

are  being  exchanged,  a  scheme  similar  to  the  exchange  jump  on  the 
6600  could  be  used  to  make  the  time  nT^  instead  of  2m  as  might  be 
supposed . 

The  exchange  time,  must  be  less  then,  or  equal  to,  the 

time  for  the  fastest  of  the  boxes  (V,  F,  or  P)  to  process  n  words. 
The  time  for  each  of  the  boxes  to  process  one  mesh  point  (one  word) 


i 
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can  be  eatimatod  by  dividing  the  total  process  time  by  N  .  The 
ainianui  time  per  word  is  that  of  the  P  box  and  ia  13t^.  Therefore, 


nT  £  (13t  )n  or 
1  o 


t  <  13t 

o 


ina  cost  of  the  aeiaory  is 


3 

C  «  c  (224N  ),  or 
M  o 


*  c  (64n)  +  c  (32o), 
o  1 


»i  tv* 

3 

2n  ♦  ra  »  7N 

c  «  cost  per  bit  of  a  fast  memory  (*  £  ^t  ) 

o  o  o 

c  *  cost  per  bit  of  r.  slow  memory  (t  £  13t  ). 

\  1  o 

Note  that  the  memory  cycle  times  can  be  longer  if  the  memory  is 

broken  up  into  modules.  If  the  primary  memory  has  four  modules,  then 

^  t  ,  but,  if  thore  are  eight  modulus,  then  wo  cannot  let  t  -  2t 
co  o  o 

because  it  world  still  take  2t  to  get  two  words,  and  the  processor 

o 

*  >uld  have  to  wait.  However,  there  is  a  bigger  advantage  in  using 
K>dulos  in  the  secondary  memory.  If  there  are  l  modules  in  the 
secMdary  memory,  then  (n/t)  words  are  transmitted  from  each  module 
with  cycle  time  t  )  and 


(  M  t  £  (13t  )n 

W  1 


t  £  13£t 
1  o 
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This  assumes  the  secondary  memory  is  a  random  access  memory.  II  the 

secondary  memory  is  not  a  random  access  memory  but  a  sequential  access 

memory  (drum,  disk,  or  tapo  drive),  then  tho  total  transfer  time  is 

the  access  timo,  t  ,  plus  the  actual  data  transfer  time  and,  to 
a 

obviate  the  processor  waiting, 


(O'. 


t  +  i  —  i  t  £  (13t  )n 
a  11/  1  o 


or 


(i)  ■  m 


£  13 1 
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Appendix  B 


A  POSSIBLE  APPROACH  TO  TURBULENCE  PROBLEMS  ON  ILLIAC  IV 


As  indicated  in  Section  V,  one  of  the  main  problemo  of  doing  fine- 

scale  turbulence  problems  on  the  I  Iliac  is  that  such  a  large  number  of 

variables  is  needed.  There  we  usually  considered  64-bit  words,  but  even 
with  32-bit  words  the  problem  is  essentially  the  same:  Many  more  numbers 
are  needed  than  can  be  held  in  the  array  memory •  Accordingly ,  numbers 
have  to  be  passed  in  and  out  of  the  I  Iliac  disk  several  times.  The 

problem  is  to  put  numbers  to  be  stored  after  an  intermediate  calcula¬ 

tion  onto  the  disk  so  that  those  numbers  which  are  to  be  used  together 
in  the  next  sequence  of  calculations  are  stored  on  the  same  or  adjacent 
pages.  The  problem  can  be  illustrated  by  considering  a  calculation 
along  the  lines  of  Method  1  of  Section  IV.  Imagine  labeling  memory 
positions  by  grid  position  so  that  neighboring  memory  units  have 
information  from  neighboring  grid  points.  At  a  given  time  step  we 
first  compute  velocity  derivatives.  Since  these  are  computed  using 
nearby  values  of  velocities,  we  can  compute  these  so  that  the  values 
of  these  at  neighboring  points  can  be  read  onto  the  disk  on  the 
appropriate  pages.  At  the  next  step  we  are  to  solve  the  Poisson 
equation  for  the  pressure#  Imagine  this  is  a  pipe  problem.  We  first 
might  Fourier  transform  in  9  and  z.  These  numbers  then  must  go  to 
the  disk.  If  these  numbers  are  denoted  by  s  (r),  what  ccmes  out  at 
a  given  time  are  all  those  for  a  fixed  r.  However,  for  the  calculation 
of  a  transform  in  r  at  the  next  stage,  we  want  all  s^(r)  for  fixed  (i,j) 
and  different  r  all  on  the  same  page.  How  do  we  do  this?  Clearly  this 
approach  is  not  acceptable,  but  this  may  not  preclude  doing  the  calculation 
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by  another  method  *:hat  does  admit  efficient  access  and  egress  from 
the  disk.  This  may  not  be  impossible# 


To  see  this,  we  -sketch  a  possible  method  of  calculation  which, 

while  certainly  not  optimal,  may  pex'haps  work.  Imagine  that  to  solve 

the  Poisson  equation  we  decompose  P  into  P  and  P_  ,  where  P,  is  to  be 

S  A  A 

a  solution  oi  the  Laplace  equation  and  P  a  particular  solution  of  the 

s 

Poisson  equation  corresponding  to  a  correct  condition  at  one  part  of 

the  boundary  and  an  assumed  condition  at  the  next  layer  of  points 

within  the  grid.  P  is  found  by  taking  transforms  in  0  and  z  for 

s 

several  planes  starting  at  a  boundary  surface.  Then  the  transform 
s 

P^  (r)  is  found  essentially  by  solving  a  second  order  difference 

equation  in  r  with  prescribed  initial  conditions  and  slope.  After 

s 

this  is  found,  the  inverse  P  (r)  is  found  by  inveiting  the  transforms. 

Now  the  pressures  at  neighboring  points  are  being  produced  together 

s 

and  can  be  efficiently  stored  in  the  disk.  When  P  has  been  calculated 
for  all  r  we  still  must  satisfy  the  correct  boundary  conditions.  This 
is  done  by  finding  a  solution  to  the  Laplace  equation  with  boundary 
values  which  are  the  difference  of  the  correct  ones  and  those  implied 

by  our  particular  solution  P  . 

s 
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